ED50 calculation
Marzonie et al.
1/08/2022
1. Experimental metadata
library(nlme)
library(aomisc)## Warning: package 'carData' was built under R version 4.0.5library(emmeans)
library(ggpubr)
library(tidyverse)## Warning: package 'tidyr' was built under R version 4.0.5## Warning: package 'readr' was built under R version 4.0.5load(file = "1_data/psii_metadata.RData")
psii = psii %>% mutate(Species = factor(Species, levels = c("A. cf humilis", "P. meandrina", "P. verrucosa")))
head(psii)## Region Reef Site Date Species Depth
## 1 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis 10.967742
## 2 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis 10.967742
## 3 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis 8.387097
## 4 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis 8.387097
## 5 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis 8.387097
## 6 Northern CSMP Bougainville Bgnv_1 2020-03-10 A. cf humilis 8.387097
## Bleaching catBleaching Treatment Tank Rack Clip Vial max.yield
## 1 2 Poor T0 2-2 24 530 A hum 275 0.471
## 2 2 Poor T3 3-4 41 233 A hum 275 0.470
## 3 3 Fair T0 2-2 24 1002 A hum 276 0.604
## 4 3 Fair T3 3-4 41 265 A hum 276 0.624
## 5 3 Fair T6 4-3 45 304 A hum 276 0.429
## 6 3 Fair T9 1-4 5 992 A hum 276 0.235
## mean.yield median.yield sd.yield TargetTemp Temperature sd.Temperature
## 1 0.4346667 0.430 0.034239354 29.99 30.57 0.05
## 2 0.4700000 0.470 0.000000000 31.81 32.06 0.24
## 3 0.5773333 0.587 0.032593455 29.99 30.57 0.05
## 4 0.6203333 0.620 0.003511885 31.81 32.06 0.24
## 5 0.4290000 0.429 NA 34.81 34.97 0.36
## 6 0.1950000 0.219 0.056000000 37.81 36.33 0.90
## MMM MMM_dev DHW duration int_mean int_max int_cum
## 1 28.96 1.61 10 32 1.386897 2.007517 12.48207
## 2 28.96 3.10 10 32 1.386897 2.007517 12.48207
## 3 28.96 1.61 10 32 1.386897 2.007517 12.48207
## 4 28.96 3.10 10 32 1.386897 2.007517 12.48207
## 5 28.96 6.01 10 32 1.386897 2.007517 12.48207
## 6 28.96 7.87 10 32 1.386897 2.007517 12.48207We’ve combined data from the sampling of coral nubbins, experimental data and results, with long-term environmental data into a single data frame for analysis. The data include the following information:
Geographic metadata Region: Bioregions reflective of the biogeoraphy of coral and fish diversity in the Coral Sea Marine Park Reef: Distinct coral reef groups Site: Sampling site
Sampling metadata Date: Date of collection Species: Species name. Pocillopora confirmed using RFLP Depth: Collection depth Bleaching: Bleaching category 1 (bleached) to 6 (healthy) catBleaching: Bleaching category grouping (Poor: 1&2, Fair: 3&4, Good:5&6)
Experimental metadata Treatment: Temperature treatment Tank: Tank within treatment jackets, 3 tanks per jacket Rack: Rack within Tank Clip: Clip on rack Vial: Unique identification number for each coral genet
Experimental measurements max.yield: maximum measurement of PSII activity from 3 measures mean.yield: mean measurement of PSII activity from 3 measures median.yield: median measurement of PSII activity from 3 measures sd.yield: standard deviation of PSII activity from 3 measures
Treatment data TargetTemp: Target treatment temperature Temperature: Average measured treatment temperature sd.Temperature: Standard deviation of measured treatment temperature
NOAA Coral Reef Watch metadata MMM: Local Maximum Monthly Mean sea surface temperature at the time of sampling MMM_dev: Deviation from MM at the time of sampling DHW: Degree Heating Weeks at the time of sampling
RmarineHeatWaves metadata (based on NOAA CRW SST data) duration: Heatwave duration (days) int_mean: Mean intensity int_max: Max intensity int_cum: Cumulative intensity
2. Dose Response curves
In order to run a log logistic regression on the experimental data and identify differences between species we need to make sure all species have a common set of reefs and enough samples from each reef and bleaching categories to explore the effects of any interactions. So we’ll create two new data frames, one for model exploration/selection (psii.1) and another for the comparison between species (psii.2). These exclude some reefs where sample sizes were low.
#How many samples from each Species/Reef. This will include all fragments (sum of all fragments, NOT individual genotypes)
psii %>% group_by(Species, Reef) %>%
dplyr::summarise(n = dplyr::n()) %>%
spread(key = Species, value = n)## `summarise()` has grouped output by 'Species'. You can override using the
## `.groups` argument.## # A tibble: 9 × 4
## Reef `A. cf humilis` `P. meandrina` `P. verrucosa`
## <fct> <int> <int> <int>
## 1 Bougainville 73 36 59
## 2 Chilcott 70 39 11
## 3 Flinders 112 20 97
## 4 Frederick 62 48 6
## 5 Herald 62 16 36
## 6 Lihou 117 23 71
## 7 Moore 67 63 37
## 8 Saumarez NA 85 NA
## 9 Wreck 96 52 24#Creating a subset of the data for model selection
psii.1 = psii %>% filter(! Reef %in% c("Saumarez", "Frederick", "Chilcott", "Wreck")) %>%
mutate(Reef = factor(Reef)) %>%
mutate(Species = factor(Species)) %>%
droplevels() %>%
as.data.frame()
#Creating a subset of the data for model selection
psii.2 = psii %>% filter(! Reef %in% c("Saumarez", "Frederick")) %>%
mutate(Reef = factor(Reef)) %>%
mutate(Species = factor(Species)) %>%
droplevels() %>%
as.data.frame()#Table 1. Number of individual corals collected from each reef AND used in the model.
#count of how many individuals were collected per reef
table1 = psii %>% group_by(Species, Reef, Vial) %>%
dplyr::summarise() %>%
group_by(Species, Reef) %>%
dplyr::summarise(ln = n()) %>%
spread(Species, ln) ## `summarise()` has grouped output by 'Species', 'Reef'. You can override using
## the `.groups` argument.
## `summarise()` has grouped output by 'Species'. You can override using the
## `.groups` argument.table1## # A tibble: 9 × 4
## Reef `A. cf humilis` `P. meandrina` `P. verrucosa`
## <fct> <int> <int> <int>
## 1 Bougainville 23 10 17
## 2 Chilcott 18 10 3
## 3 Flinders 30 5 26
## 4 Frederick 17 13 2
## 5 Herald 18 5 11
## 6 Lihou 32 6 18
## 7 Moore 20 17 10
## 8 Saumarez NA 22 NA
## 9 Wreck 24 13 6psii %>% group_by(Species, Reef, Treatment) %>%
summarise(mean = round(mean(median.yield),2)) %>%
spread(key = Treatment, value = mean)## `summarise()` has grouped output by 'Species', 'Reef'. You can override using
## the `.groups` argument.## # A tibble: 25 × 6
## # Groups: Species, Reef [25]
## Species Reef T0 T3 T6 T9
## <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 A. cf humilis Bougainville 0.51 0.55 0.45 0.1
## 2 A. cf humilis Chilcott 0.62 0.63 0.47 0.07
## 3 A. cf humilis Flinders 0.61 0.61 0.44 0.08
## 4 A. cf humilis Frederick 0.62 0.62 0.51 0.08
## 5 A. cf humilis Herald 0.6 0.59 0.35 0.04
## 6 A. cf humilis Lihou 0.61 0.61 0.49 0.03
## 7 A. cf humilis Moore 0.55 0.57 0.41 0.07
## 8 A. cf humilis Wreck 0.65 0.65 0.61 0.18
## 9 P. meandrina Bougainville 0.61 0.61 0.58 0.1
## 10 P. meandrina Chilcott 0.64 0.65 0.56 0.05
## # … with 15 more rows
## # ℹ Use `print(n = ...)` to see more rowspsii %>% group_by(Species, Reef, Treatment) %>%
summarise(mean = round(sd(median.yield),2)) %>%
spread(key = Treatment, value = mean)## `summarise()` has grouped output by 'Species', 'Reef'. You can override using
## the `.groups` argument.## # A tibble: 25 × 6
## # Groups: Species, Reef [25]
## Species Reef T0 T3 T6 T9
## <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 A. cf humilis Bougainville 0.08 0.06 0.07 0.07
## 2 A. cf humilis Chilcott 0.03 0.02 0.08 0.05
## 3 A. cf humilis Flinders 0.04 0.03 0.13 0.1
## 4 A. cf humilis Frederick 0.02 0.03 0.18 0.08
## 5 A. cf humilis Herald 0.03 0.04 0.09 0.04
## 6 A. cf humilis Lihou 0.03 0.03 0.13 0.02
## 7 A. cf humilis Moore 0.06 0.05 0.13 0.08
## 8 A. cf humilis Wreck 0.02 0.02 0.05 0.11
## 9 P. meandrina Bougainville 0.08 0.07 0.07 0.08
## 10 P. meandrina Chilcott 0.02 0.04 0.07 0.04
## # … with 15 more rows
## # ℹ Use `print(n = ...)` to see more rowsResults 3.2 Differences in Fv/Fm observed in different treatments
#reporting differences among treatments in results section 3.2
psii.2 %>%
dplyr::group_by(Treatment) %>%
summarise(med.treat = mean(median.yield), sd.treat = sd(median.yield), median.treat = median(median.yield))## # A tibble: 4 × 4
## Treatment med.treat sd.treat median.treat
## <fct> <dbl> <dbl> <dbl>
## 1 T0 0.610 0.0567 0.624
## 2 T3 0.619 0.0486 0.628
## 3 T6 0.546 0.114 0.588
## 4 T9 0.0855 0.0899 0.05Assuming a fixed lower asymptote representing complete loss of photosynthetic yield, a three parameter log-logistic model can be assumed, estimating the upper asymptote, the location of the inflection point (PSII-50), and steepness of the curve whereby:
b is the slope around the inflection point d is the upper asymptote, e is the inflection point (PSII-50)
2.1 Model selection : SPECIES
What parameters to estimate?
Our primary interest is in the inflection point (e or ED50) but other parameter estimates may be important and may vary among species. We anticipate most of the variation in a species response to temperature treatment at this scale will occur at the Reef level so we included Reef as a random factor to explore the importance of parameter estimates
#Naive model to identify starting points
# Model Parameters: b = shape parameter (i.e steepness); d=upper limit; e=ED50
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.1,
curveid = Species,
pmodels = c( ~ 1, ~ 1, ~ Species), bcVal = 0.5)
species.1 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
random = e ~ 1|Reef,
fixed = list(b ~ 1, d ~ 1, e ~ Species),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.1,
curveid = Species,
pmodels = c( ~ 1, ~ Species, ~ Species), bcVal = 0.5)
species.2 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
random = e ~ 1|Reef,
fixed = list(b ~ 1, d ~ Species, e ~ Species),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.1,
curveid = Species,
pmodels = c( ~ Species, ~ 1, ~ Species), bcVal = 0.5)
species.3 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
random = e ~ 1|Reef,
fixed = list(b ~ Species, d ~ 1, e ~ Species),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.1,
curveid = Species,
pmodels = c( ~ Species, ~ Species, ~ Species), bcVal = 0.5)
species.4 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
random = e ~ 1|Reef,
fixed = list(b ~ Species, d ~ Species, e ~ Species),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
AIC(species.1,species.2,species.3, species.4)## df AIC
## species.1 7 -2131.489
## species.2 9 -2179.997
## species.3 9 -2171.099
## species.4 11 -2203.347The more complex model that estimates all 3 parameters (species.4) is the best model. However, it may be necessary to reduce the complexity of the model when accounting for other random factors in the model, in which case estimating d and e provide the best result (species.2). For example, when accounting for Depth and Tank effects, the models are too complex to also estimate b.
Accounting for Depth effects
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.1,
curveid = Species,
pmodels = c( ~ 1, ~ Species, ~ Species), bcVal = 0.5)
species.5 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
random = e ~ 1|Reef/Depth,
fixed = list(b ~ 1, d ~ Species, e ~ Species),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
AIC(species.2,species.5) ## df AIC
## species.2 9 -2179.997
## species.5 10 -2177.997anova(species.2, species.5)## Model df AIC BIC logLik Test L.Ratio p-value
## species.2 1 9 -2179.997 -2136.886 1098.998
## species.5 2 10 -2177.997 -2130.096 1098.998 1 vs 2 8.668175e-06 0.9977Depth does not need to be accounted for.
Accounting for Tank effects
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.1,
curveid = Species,
pmodels = c( ~ 1, ~ Species, ~ Species), bcVal = 0.5)
species.6 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
random = e ~ 1|Reef/Tank,
fixed = list(b ~ 1, d ~ Species, e ~ Species),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
AIC(species.2,species.6)## df AIC
## species.2 9 -2179.997
## species.6 10 -2177.081anova(species.2, species.6)## Model df AIC BIC logLik Test L.Ratio p-value
## species.2 1 9 -2179.997 -2136.886 1098.998
## species.6 2 10 -2177.081 -2129.180 1098.541 1 vs 2 0.9153679 0.3387Tank does not need to be accounted for.
Accounting for bleaching
Introduce an interaction between Species and Bleaching categories. We explore the interaction only for d and e.
#effect of bleached samples on e
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.1,
curveid = Species:catBleaching,
pmodels = c( ~ 1, ~ Species, ~ Species*catBleaching), bcVal = 0.5)
species.7 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
random = e ~ 1|Reef,
fixed = list(b ~ 1, d ~ Species, e ~ Species*catBleaching),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
#effect of bleached samples on d
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.1,
curveid = Species:catBleaching,
pmodels = c( ~ 1, ~ Species*catBleaching, ~ Species), bcVal = 0.5)
species.8 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
random = e ~ 1|Reef,
fixed = list(b ~ 1, d ~ Species*catBleaching, e ~ Species),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
#effect of bleached samples on d and e
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.1,
curveid = Species:catBleaching,
pmodels = c( ~ 1, ~ Species*catBleaching, ~ Species*catBleaching), bcVal = 0.5)
species.9 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
random = e ~ 1|Reef,
fixed = list(b ~ 1, d ~ Species*catBleaching, e ~ Species*catBleaching),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
#check random factors on d and e
species.9.1 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.1,
random = d + e ~ 1|Reef,
fixed = list(b ~ 1, d ~ Species*catBleaching, e ~ Species*catBleaching),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.1, : Iteration 1, LME step: nlminb() did not converge (code = 1). PORT
## message: function evaluation limit reached without convergence (9)AIC(species.2, species.7, species.8, species.9, species.9.1)## df AIC
## species.2 9 -2179.997
## species.7 15 -2181.960
## species.8 15 -2230.468
## species.9 21 -2244.262
## species.9.1 23 -2277.888anova(species.2,species.9.1)## Model df AIC BIC logLik Test L.Ratio p-value
## species.2 1 9 -2179.997 -2136.886 1098.998
## species.9.1 2 23 -2277.888 -2167.716 1161.944 1 vs 2 125.8915 <.0001anova(species.9.1) # interaction is significant## numDF denDF F-value p-value
## b 1 866 875.517 <.0001
## d.(Intercept) 1 866 315.787 <.0001
## d.Species 2 866 26.764 <.0001
## d.catBleaching 2 866 117.710 <.0001
## d.Species:catBleaching 4 866 9.038 <.0001
## e.(Intercept) 1 866 24674.744 <.0001
## e.Species 2 866 89.599 <.0001
## e.catBleaching 2 866 10.094 <.0001
## e.Species:catBleaching 4 866 2.932 0.02The interaction between bleaching categories and species is important and improves the model so we need to control for the level of bleaching of sampled nubbins. The effect on e and d is small but should not be ignored. If we look at the contrast below, the species differences change depending on the bleaching category.
#The effect of including the interaction between species and catBleaching on d and e is small
emmeans(species.9.1,pairwise~catBleaching|Species, param="d", level = 0.95)$emmeans %>% plot(comparisons = TRUE)
emmeans(species.9.1,pairwise~catBleaching|Species, param="e", level = 0.95)$emmeans %>% plot(comparisons = TRUE)
#The contrast show how bleaching can affect the species comparisons
emmeans(species.9.1,pairwise~catBleaching|Species, param="d", level = 0.95)$contrast %>% plot()
emmeans(species.9.1,pairwise~Species|catBleaching, param="d", level = 0.95)$contrast %>% plot()
emmeans(species.9.1,pairwise~catBleaching|Species, param="e", level = 0.95)$contrast %>% plot()
emmeans(species.9.1,pairwise~Species|catBleaching, param="e", level = 0.95)$contrast %>% plot()
Note: The blue bars are confidence intervals for the EMMs, and the red
arrows are for the comparisons among them. If an arrow from one mean
overlaps an arrow from another group, the difference is not
“significant,” based on the adjust setting (which defaults to “tukey”)
and the value of alpha (which defaults to 0.05). (https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html).
2.2 Model predictions: Species*Bleaching
#including larger subset of reefs (psii.2)
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.2,
curveid = Species:catBleaching,
pmodels = c( ~ 1, ~ Species*catBleaching, ~ Species*catBleaching), bcVal = 0.5)
species.9.2 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.2,
random = d + e ~ 1|Reef,
fixed = list(b ~ 1, d ~ Species*catBleaching, e ~ Species*catBleaching),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.2, : Iteration 1, LME step: nlminb() did not converge (code = 1). PORT
## message: function evaluation limit reached without convergence (9)#Chilcott and Wreck have a strong effect on parameter estimates of e and d for each species particularly for A. cf humilis
emmeans(species.9.1, ~Species, param = 'e')## NOTE: Results may be misleading due to involvement in interactions## Species emmean SE df lower.CL upper.CL
## A. cf humilis 6.76 0.0607 866 6.64 6.87
## P. meandrina 7.33 0.0814 866 7.17 7.49
## P. verrucosa 7.73 0.0663 866 7.60 7.86
##
## Results are averaged over the levels of: catBleaching
## Confidence level used: 0.95emmeans(species.9.2, ~Species, param = 'e')## NOTE: Results may be misleading due to involvement in interactions## Species emmean SE df lower.CL upper.CL
## A. cf humilis 7.10 0.140 1156 6.82 7.37
## P. meandrina 7.36 0.143 1156 7.08 7.64
## P. verrucosa 7.73 0.143 1156 7.45 8.01
##
## Results are averaged over the levels of: catBleaching
## Confidence level used: 0.95emmeans(species.9.1, ~Species, param = 'd')## NOTE: Results may be misleading due to involvement in interactions## Species emmean SE df lower.CL upper.CL
## A. cf humilis 0.578 0.00814 866 0.562 0.594
## P. meandrina 0.624 0.00970 866 0.605 0.644
## P. verrucosa 0.618 0.00845 866 0.602 0.635
##
## Results are averaged over the levels of: catBleaching
## Confidence level used: 0.95emmeans(species.9.2, ~Species, param = 'd')## NOTE: Results may be misleading due to involvement in interactions## Species emmean SE df lower.CL upper.CL
## A. cf humilis 0.585 0.00614 1156 0.573 0.597
## P. meandrina 0.632 0.00747 1156 0.618 0.647
## P. verrucosa 0.628 0.00692 1156 0.615 0.642
##
## Results are averaged over the levels of: catBleaching
## Confidence level used: 0.95#what about the slope b
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.2,
curveid = Species:catBleaching,
pmodels = c( ~ Species*catBleaching, ~ Species*catBleaching, ~ Species*catBleaching), bcVal = 0.5)
species.9.3 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.2,
random = b + d + e ~ 1|Reef,
fixed = list(b ~ Species*catBleaching, d ~ Species*catBleaching, e ~ Species*catBleaching),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.2, : Iteration 1, LME step: nlminb() did not converge (code = 1). PORT
## message: function evaluation limit reached without convergence (9)## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.2, : Singular precision matrix in level -1, block 1
## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.2, : Singular precision matrix in level -1, block 1
## Warning in nlme.formula(median.yield ~ NLS.L3(MMM_dev, b, d, e), data =
## psii.2, : Singular precision matrix in level -1, block 1modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.2,
curveid = Species:catBleaching,
pmodels = c( ~ Species, ~ Species*catBleaching, ~ Species*catBleaching), bcVal = 0.5)
#Does not run
#species.9.4 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.2,
# random = b + d + e ~ 1|Reef,
# fixed = list(b ~ Species, d ~ Species*catBleaching, e ~ Species*catBleaching),
# start = coef(modNaive.Species),
# control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
modNaive.Species <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.2,
curveid = Species:catBleaching,
pmodels = c( ~ catBleaching, ~ Species*catBleaching, ~ Species*catBleaching), bcVal = 0.5)
species.9.5 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.2,
random = b + d + e ~ 1|Reef,
fixed = list(b ~ catBleaching, d ~ Species*catBleaching, e ~ Species*catBleaching),
start = coef(modNaive.Species),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
AIC(species.9.2, species.9.3, species.9.5)## df AIC
## species.9.2 23 -2982.373
## species.9.3 34 -3109.717
## species.9.5 28 -3034.841anova(species.9.2, species.9.3)## Model df AIC BIC logLik Test L.Ratio p-value
## species.9.2 1 23 -2982.373 -2865.669 1514.187
## species.9.3 2 34 -3109.718 -2937.198 1588.859 1 vs 2 149.3442 <.0001Species also have different slopes and so it is important that we consider this in our model. Significant improvements to model fit with species.9.3
Estimated marginal mean & pairwise contrasts
parameter e
my.model1 = species.9.3
load("3_outputs/nlme_All species.RData")
my.data1 = psii.2
emmeans(my.model1, ~Species, param = 'e', level = 0.95) %>% plot(comparisons = TRUE)## NOTE: Results may be misleading due to involvement in interactions
emmeans(my.model1, pairwise~Species, param = 'e', level = 0.95)$contrast %>% plot()## NOTE: Results may be misleading due to involvement in interactions
The differences in the inflection point (e, ED50) between species are
significant.
emmeans(my.model1, pairwise~Species, param = 'e', level = 0.95)$contrast #ED50 contrasts ## NOTE: Results may be misleading due to involvement in interactions## contrast estimate SE df t.ratio p.value
## A. cf humilis - P. meandrina -0.376 0.0765 1148 -4.914 <.0001
## A. cf humilis - P. verrucosa -0.696 0.0767 1148 -9.067 <.0001
## P. meandrina - P. verrucosa -0.320 0.0857 1148 -3.733 0.0006
##
## Results are averaged over the levels of: catBleaching
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(my.model1, ~Species, param = 'e', level = 0.95) ## NOTE: Results may be misleading due to involvement in interactions## Species emmean SE df lower.CL upper.CL
## A. cf humilis 7.05 0.154 1148 6.75 7.35
## P. meandrina 7.42 0.159 1148 7.11 7.74
## P. verrucosa 7.74 0.160 1148 7.43 8.06
##
## Results are averaged over the levels of: catBleaching
## Confidence level used: 0.95parameter d
emmeans(my.model1, ~Species, param = 'd', level = 0.95) %>% plot(comparisons = TRUE)## NOTE: Results may be misleading due to involvement in interactions
emmeans(my.model1, pairwise~Species, param = 'd', level = 0.95)$contrast %>% plot()## NOTE: Results may be misleading due to involvement in interactions
A. cf humilis has lower upper asymptote than either P. meandrina or P.
verrucosa. There is no observable difference between Pocilloporids.
emmeans(my.model1, pairwise~Species, param = 'd', level = 0.95)$contrast #ED50 contrasts ## NOTE: Results may be misleading due to involvement in interactions## contrast estimate SE df t.ratio p.value
## A. cf humilis - P. meandrina -0.03078 0.00696 1148 -4.421 <.0001
## A. cf humilis - P. verrucosa -0.03531 0.00657 1148 -5.378 <.0001
## P. meandrina - P. verrucosa -0.00453 0.00781 1148 -0.581 0.8306
##
## Results are averaged over the levels of: catBleaching
## P value adjustment: tukey method for comparing a family of 3 estimatesparameter b
emmeans(my.model1, ~Species, param = 'b', level = 0.95) %>% plot(comparisons = TRUE)## NOTE: Results may be misleading due to involvement in interactions
emmeans(my.model1, pairwise~Species, param = 'b', level = 0.95)$contrast %>% plot()## NOTE: Results may be misleading due to involvement in interactions
A. cf humilis has slacker slope than either P. meandrina or P.
verrucosa. The condition of A. cf humilis may have dropped at lower
temps. There is no observable difference between Pocilloporids.
emmeans(my.model1, pairwise~Species, param = 'b', level = 0.95)$contrast #ED50 contrasts ## NOTE: Results may be misleading due to involvement in interactions## contrast estimate SE df t.ratio p.value
## A. cf humilis - P. meandrina 0.853 0.156 1148 5.455 <.0001
## A. cf humilis - P. verrucosa 0.663 0.178 1148 3.722 0.0006
## P. meandrina - P. verrucosa -0.190 0.218 1148 -0.875 0.6565
##
## Results are averaged over the levels of: catBleaching
## P value adjustment: tukey method for comparing a family of 3 estimatescatBleaching interaction
#Including bleaching
emmeans(my.model1, ~catBleaching|Species, param = 'e', level = 0.95) %>% plot()
emmeans(my.model1, ~catBleaching|Species, param = 'd', level = 0.95) %>% plot()
emmeans(my.model1, ~catBleaching|Species, param = 'b', level = 0.95) %>% plot()
catBleaching has its strongest effect on d within species, whereby
colonies that started in a poorer condition have lower values of median
PSII yield.
SOM Tables S4 and S4: Species Pairwise Contrasts and Emmeans
emmeans(my.model1, pairwise~Species, param = 'e', level = 0.95)$contrast #ED50 contrasts ## NOTE: Results may be misleading due to involvement in interactions## contrast estimate SE df t.ratio p.value
## A. cf humilis - P. meandrina -0.376 0.0765 1148 -4.914 <.0001
## A. cf humilis - P. verrucosa -0.696 0.0767 1148 -9.067 <.0001
## P. meandrina - P. verrucosa -0.320 0.0857 1148 -3.733 0.0006
##
## Results are averaged over the levels of: catBleaching
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(my.model1, pairwise~Species, param = 'd', level = 0.95)$contrast #ED50 contrasts ## NOTE: Results may be misleading due to involvement in interactions## contrast estimate SE df t.ratio p.value
## A. cf humilis - P. meandrina -0.03078 0.00696 1148 -4.421 <.0001
## A. cf humilis - P. verrucosa -0.03531 0.00657 1148 -5.378 <.0001
## P. meandrina - P. verrucosa -0.00453 0.00781 1148 -0.581 0.8306
##
## Results are averaged over the levels of: catBleaching
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(my.model1, pairwise~Species, param = 'b', level = 0.95)$contrast #ED50 contrasts ## NOTE: Results may be misleading due to involvement in interactions## contrast estimate SE df t.ratio p.value
## A. cf humilis - P. meandrina 0.853 0.156 1148 5.455 <.0001
## A. cf humilis - P. verrucosa 0.663 0.178 1148 3.722 0.0006
## P. meandrina - P. verrucosa -0.190 0.218 1148 -0.875 0.6565
##
## Results are averaged over the levels of: catBleaching
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans(my.model1, ~Species, param = 'e', level = 0.95) ## NOTE: Results may be misleading due to involvement in interactions## Species emmean SE df lower.CL upper.CL
## A. cf humilis 7.05 0.154 1148 6.75 7.35
## P. meandrina 7.42 0.159 1148 7.11 7.74
## P. verrucosa 7.74 0.160 1148 7.43 8.06
##
## Results are averaged over the levels of: catBleaching
## Confidence level used: 0.95Model check
resid(my.model1) %>% plot()
VarCorr(my.model1)## Reef = pdLogChol(list(b ~ 1,d ~ 1,e ~ 1))
## Variance StdDev Corr
## b.(Intercept) 2.979937e-20 1.726249e-10 b.(In) d.(In)
## d.(Intercept) 4.497247e-04 2.120671e-02 -0.662
## e.(Intercept) 1.514108e-01 3.891154e-01 -0.935 0.374
## Residual 3.829744e-03 6.188493e-02ranef(my.model1)## b.(Intercept) d.(Intercept) e.(Intercept)
## Bougainville 1.539353e-10 -0.036360850 -0.130969244
## Chilcott 3.784724e-11 0.022342256 -0.293119010
## Flinders -7.389782e-11 -0.003692153 0.239259916
## Herald 1.733321e-10 -0.001971101 -0.472525825
## Lihou -1.780246e-11 0.006769950 -0.006246521
## Moore 9.070098e-11 -0.015226812 -0.128861706
## Wreck -3.641153e-10 0.028138711 0.792462390#coef(my.model1)
#vcov(my.model1)2.3 Figure: Species comparison
##Create predictions data frame
predframe.spp = with(my.data1,
expand.grid(Species = levels(Species),
catBleaching = levels(catBleaching),
MMM_dev=seq(min(MMM_dev), max(MMM_dev) + 1, 0.1)))
#Adding confidence intervals
prdc.spp = as.data.frame(nlraa::predict_nlme(my.model1, newdata = predframe.spp,
interval = "confidence", level = 0.95))
prdc.spp$Species = predframe.spp$Species
prdc.spp$MMM_dev = predframe.spp$MMM_dev
prdc.spp$catBleaching = predframe.spp$catBleaching
#Averaging predictions across catBleaching
prdc.average = prdc.spp %>% dplyr::select(-catBleaching) %>%
group_by(Species, MMM_dev) %>% summarise_all(.funs=mean)
# Predict ED50 for each species
ED50spp = emmeans(my.model1, ~Species|catBleaching, param= "e") %>% as_tibble()
predframe.ED50 = ED50spp %>% dplyr::select(Species, catBleaching, emmean) %>%
dplyr::rename(MMM_dev = emmean)
prdc.ED50 = as.data.frame(nlraa::predict_nlme(my.model1, newdata = predframe.ED50,
interval = "confidence", level = 0.95))
prdc.ED50 = bind_cols(predframe.ED50, prdc.ED50)
prdc.ED50.average = prdc.ED50 %>% dplyr::select(-catBleaching) %>%
group_by(Species) %>% summarise_all(.funs=mean)Fig 2A: PSII comparison by species
library(wesanderson)
col = wes_palette("GrandBudapest1", 4, type = "continuous")
Species.plot = ggplot() +
geom_segment(data = prdc.ED50.average, aes(x = MMM_dev, xend = MMM_dev, y = -Inf,
yend = Estimate , color = Species),
linetype = 1, size = .3) +
geom_point(data=psii.1, aes(x=MMM_dev, y=median.yield, col=Species), alpha=0.1, size = 1) +
geom_line(data = prdc.average, aes(colour=Species, y=Estimate, x=MMM_dev), size = .3) +
geom_ribbon(data=prdc.average, aes(x = MMM_dev, ymin = Q2.5, ymax = Q97.5, fill = Species),
alpha = 0.5, show.legend = FALSE) +
geom_point(data = prdc.ED50.average, aes(x = MMM_dev, y = Estimate, fill = Species),
size = 1.5, stroke = .2, shape = 21, col = "white") +
scale_x_continuous(breaks = c(3,6,9)) +
coord_cartesian(xlim = c(2,9.6), ylim = c(0,.75)) +
scale_fill_manual(values=col) +
scale_colour_manual(values=col) +
labs(x="Temp. above local MMM (°C)",
y="F[v]/F[m]") +
theme_classic() +
theme(axis.line = element_blank(),
panel.border = element_rect(size = .5, fill = "transparent"),
legend.title = element_blank(),
legend.background = element_rect(fill = "transparent"),
#legend.position = "none",
legend.text = element_text(face = "italic", size = 7, family = "Helvetica"),
legend.key.size = unit(2, units = "mm"),
legend.margin = margin(0,0,0,0, unit = "mm"),
axis.text.y = element_text(size = 8, family = "Helvetica"),
axis.title = element_text(size = 9, family = "Helvetica"))
Species.plot
ggsave("2_figures/SpeciesPlotonly.png", width = 140, height = 70, units = "mm")Fig SX: PSII Species X catBleaching
SpeciesBl.plot = ggplot() +
geom_segment(data = prdc.ED50, aes(x = MMM_dev, xend = MMM_dev, y = -Inf, yend = Estimate , color = Species),
linetype = 1, size = .3) +
geom_point(data=psii.1, aes(x=MMM_dev, y=median.yield, col=Species), alpha=0.1, size = 1) +
geom_line(data = prdc.spp, aes(colour=Species, y=Estimate, x=MMM_dev), size = .3) +
geom_ribbon(data=prdc.spp, aes(x = MMM_dev, ymin = Q2.5, ymax = Q97.5, fill = Species),
alpha = 0.5, show.legend = FALSE) +
geom_point(data = prdc.ED50, aes(x = MMM_dev, y = Estimate, fill = Species), size = 1.5, stroke = .2, shape = 21, col = "white") +
scale_x_continuous(breaks = c(3,6,9)) +
coord_cartesian(xlim = c(2,9.6), ylim = c(0,.75)) +
scale_fill_manual(values=col) +
scale_colour_manual(values=col) +
labs(x="Temp. above local MMM (°C)",
y="Fv/Fm") +
theme_classic() +
facet_wrap(~catBleaching) +
theme(axis.line = element_blank(),
panel.border = element_rect(size = .5, fill = "transparent"),
legend.title = element_blank(),
legend.background = element_rect(fill = "transparent"),
legend.position = "right",
legend.text = element_text(face = "italic", size = 7, family = "Helvetica"),
legend.key.size = unit(2, units = "mm"),
legend.margin = margin(0,0,0,0, unit = "mm"),
axis.text.y = element_text(size = 8, family = "Helvetica"),
axis.title = element_text(size = 9, family = "Helvetica"))
SpeciesBl.plot
ggsave("2_figures/FigS3_catBleaching_by_Species.png", width = 140, height = 70, units = "mm")2.4 Model selection : REEF
Combined 3-way interaction
Since we don’t have the same set of reefs for all species, we can’t run the three-way interaction Species:Reef:catBleaching but we can still explore what this will look like for a subset of reefs and how it affects each model parameter. However, this significantly increases the complexity of the model so we need to simplify it.
If we look at the interaction between species and catBleaching. We saw that for species, the most complex model, with an interaction between species and catBleaching on all parameter estimates was important. However, the strength of the effect depends on the species and parameter. We can see from below: - catBleaching has a weak effect on b, though can be an important source of variation between species - catBleaching has a strong effect on d, particularly for A. humilis - catBleaching has a weak effect on e, with the only significant effect on P. verrucosa
emmeans(my.model1, ~catBleaching|Species, param = 'b', level = 0.95) %>% plot(comparisons = TRUE)
emmeans(my.model1, ~catBleaching|Species, param = 'd', level = 0.95) %>% plot(comparisons = TRUE)
emmeans(my.model1, ~catBleaching|Species, param = 'e', level = 0.95) %>% plot(comparisons = TRUE)
Therefore when developing a more complex model that include reefs. We
must consider: - Species differences on b - Species and catBleahching
differences on d - Species differences on e
However, running Species on b does not work, perhaps because it is confounded by another fixed factor. So we have to fix the slope.
#testing the 3-interaction on psii.2 (a common subset of reefs)
modNaive <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = psii.2,
curveid = Species:Reef:catBleaching,
pmodels = c( ~ 1, ~ Species*catBleaching, ~ Species*Reef), bcVal = 0.5)
reef.1 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = psii.2,
random = e ~ 1|Vial,
fixed = list(b ~ 1, d ~ Species*catBleaching, e ~ Species*Reef),
start = coef(modNaive),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
emmeans(reef.1, ~Reef|Species, param = "e") %>% plot(comparisons = TRUE)## Warning: Comparison discrepancy in group "P. verrucosa", Bougainville - Wreck:
## Target overlap = 0.0101, overlap on graph = -0.0223
emmeans(reef.1, ~catBleaching|Species, param = "d") %>% plot(comparisons = TRUE)
anova(reef.1)## numDF denDF F-value p-value
## b 1 829 8603.21 <.0001
## d.(Intercept) 1 829 151066.39 <.0001
## d.Species 2 829 58.34 <.0001
## d.catBleaching 2 829 191.14 <.0001
## d.Species:catBleaching 4 829 30.84 <.0001
## e.(Intercept) 1 829 100928.49 <.0001
## e.Species 2 829 13.86 <.0001
## e.Reef 6 829 85.86 <.0001
## e.Species:Reef 12 829 11.11 <.0001##Create predictions data frame
predframe.spp = with(psii.2, expand.grid(Species = levels(Species),
Reef = levels(Reef),
catBleaching = levels(catBleaching),
MMM_dev=seq(min(MMM_dev), max(MMM_dev) + 1, 0.1)))
#Adding confidence intervals
prdc.spp = as.data.frame(nlraa::predict_nlme(reef.1, newdata = predframe.spp,
interval = "confidence", level = 0.950))
prdc.spp = bind_cols(predframe.spp, prdc.spp)
ggplot(prdc.spp, aes(x = MMM_dev, y = Estimate, col = Reef)) +
geom_line() +
facet_grid(catBleaching~Species) +
theme_light()
ggplot(prdc.spp, aes(x = MMM_dev, y = Estimate, col = catBleaching)) +
geom_line() +
facet_grid(Reef~Species) +
theme_light()
d.Species:catBleaching - The interaction is significant e.Species:Reef -
The interaction is significant
Species-specific models
Due to some species not being present for some reefs, we will run the species separately. The difficulty now is that not all bleaching categories will be present for each reef in each species. We can compare the results obtained from the 3-way interaction with the results obtained for each individuals species modelled separately.
First, we’ll have to subset the data to remove reefs with small sample sizes for each reef
#How many samples from each Species/Reef
psii %>% group_by(Species, Reef) %>%
dplyr::summarise(n = dplyr::n()) %>%
spread(key = Species, value = n)## `summarise()` has grouped output by 'Species'. You can override using the
## `.groups` argument.## # A tibble: 9 × 4
## Reef `A. cf humilis` `P. meandrina` `P. verrucosa`
## <fct> <int> <int> <int>
## 1 Bougainville 73 36 59
## 2 Chilcott 70 39 11
## 3 Flinders 112 20 97
## 4 Frederick 62 48 6
## 5 Herald 62 16 36
## 6 Lihou 117 23 71
## 7 Moore 67 63 37
## 8 Saumarez NA 85 NA
## 9 Wreck 96 52 24ahum.1 = psii %>% filter(Species == "A. cf humilis") %>%
as.data.frame() %>%
droplevels()
pmea.1 = psii %>% filter(Species == "P. meandrina") %>%
droplevels()
pver.1 = psii %>%
filter(Species == "P. verrucosa") %>%
filter(! Reef %in% c("Frederick", "Chilcott")) %>%
droplevels()What parameters to estimate?
Not all bleaching categories present at each reef, explore model selection with a subset of reefs
ahum.1 %>% group_by(Reef, catBleaching) %>%
dplyr::summarise(n = dplyr::n()) %>%
spread(key = catBleaching, value = n)## `summarise()` has grouped output by 'Reef'. You can override using the
## `.groups` argument.## # A tibble: 8 × 4
## # Groups: Reef [8]
## Reef Healthy Fair Poor
## <fct> <int> <int> <int>
## 1 Bougainville 4 34 35
## 2 Chilcott 24 40 6
## 3 Flinders 38 66 8
## 4 Frederick 25 37 NA
## 5 Herald 18 35 9
## 6 Lihou 27 72 18
## 7 Moore 11 46 10
## 8 Wreck 52 44 NAahum.2 = ahum.1 %>% filter(! Reef %in% c("Frederick", "Wreck")) %>%
as.data.frame() %>% droplevels()#model d and e?
# model with
modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = ahum.2,
curveid = Reef,
pmodels = c( ~ 1, ~ 1, ~ Reef), bcVal = 0.5)
reef.ahum.1 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = e ~ 1|Vial,
fixed = list(b ~ 1, d ~ 1, e ~ Reef),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = ahum.2,
curveid = Reef,
pmodels = c( ~ 1, ~ Reef, ~ Reef), bcVal = 0.5)
reef.ahum.2 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = e ~ 1|Vial,
fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = ahum.2,
curveid = Reef,
pmodels = c( ~ Reef, ~ Reef, ~ Reef), bcVal = 0.5)
#reef.ahum.3 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
# random = e ~ 1|Vial,
# fixed = list(b ~ Reef, d ~ Reef, e ~ Reef),
# start = coef(modNaive.ahum),
# control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
AIC(reef.ahum.1, reef.ahum.2)## df AIC
## reef.ahum.1 10 -1177.798
## reef.ahum.2 15 -1221.154Can’t run any models to estimate b, we can assume that the gradient is the same within species. However, we get marginal improvements with d and e. Best model: reef.ahum.2
Accounting for bleaching
How does the catBleaching affect e and d
modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = ahum.2,
curveid = Reef:catBleaching,
pmodels = c( ~ 1, ~ Reef*catBleaching, ~ Reef*catBleaching), bcVal = 0.5)
reef.ahum.3 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = e ~ 1|Vial,
fixed = list(b ~ 1, d ~ Reef*catBleaching, e ~ Reef*catBleaching),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = ahum.2,
curveid = Reef:catBleaching,
pmodels = c( ~ 1, ~ Reef, ~ Reef*catBleaching), bcVal = 0.5)
reef.ahum.4 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = e ~ 1|Vial,
fixed = list(b ~ 1, d ~ Reef, e ~ Reef*catBleaching),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = ahum.2,
curveid = Reef:catBleaching,
pmodels = c( ~ 1, ~ Reef*catBleaching, ~ Reef), bcVal = 0.5)
reef.ahum.5 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = e ~ 1|Vial,
fixed = list(b ~ 1, d ~ Reef*catBleaching, e ~ Reef),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
AIC(reef.ahum.2, reef.ahum.3, reef.ahum.4, reef.ahum.5)## df AIC
## reef.ahum.2 15 -1221.154
## reef.ahum.3 39 -1220.852
## reef.ahum.4 27 -1221.605
## reef.ahum.5 27 -1225.515anova(reef.ahum.2, reef.ahum.5)## Model df AIC BIC logLik Test L.Ratio p-value
## reef.ahum.2 1 15 -1221.154 -1157.905 625.5770
## reef.ahum.5 2 27 -1225.515 -1111.667 639.7575 1 vs 2 28.36111 0.0049summary(reef.ahum.5)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: median.yield ~ NLS.L3(MMM_dev, b, d, e)
## Data: ahum.2
## AIC BIC logLik
## -1225.515 -1111.667 639.7575
##
## Random effects:
## Formula: e ~ 1 | Vial
## e.(Intercept) Residual
## StdDev: 1.864301e-05 0.06748164
##
## Fixed effects: list(b ~ 1, d ~ Reef * catBleaching, e ~ Reef)
## Value Std.Error DF t-value p-value
## b -1.351935 0.05987450 336 -22.57947 0.0000
## d.(Intercept) 0.544895 0.04303904 336 12.66048 0.0000
## d.ReefChilcott 0.091942 0.04669440 336 1.96901 0.0498
## d.ReefFlinders 0.055393 0.04537504 336 1.22077 0.2230
## d.ReefHerald 0.077505 0.04804852 336 1.61305 0.1077
## d.ReefLihou 0.081456 0.04623146 336 1.76192 0.0790
## d.ReefMoore 0.068564 0.05027447 336 1.36378 0.1735
## d.catBleachingFair 0.027588 0.04510576 336 0.61162 0.5412
## d.catBleachingPoor -0.042135 0.04500108 336 -0.93630 0.3498
## d.ReefChilcott:catBleachingFair -0.035716 0.05022179 336 -0.71116 0.4775
## d.ReefFlinders:catBleachingFair -0.018777 0.04833243 336 -0.38850 0.6979
## d.ReefHerald:catBleachingFair -0.056621 0.05192923 336 -1.09036 0.2763
## d.ReefLihou:catBleachingFair -0.043119 0.04903014 336 -0.87944 0.3798
## d.ReefMoore:catBleachingFair -0.068826 0.05326301 336 -1.29219 0.1972
## d.ReefChilcott:catBleachingPoor 0.022370 0.05755472 336 0.38868 0.6978
## d.ReefFlinders:catBleachingPoor 0.043943 0.05610708 336 0.78320 0.4341
## d.ReefHerald:catBleachingPoor 0.002716 0.05683793 336 0.04778 0.9619
## d.ReefLihou:catBleachingPoor 0.028739 0.05197094 336 0.55298 0.5806
## d.ReefMoore:catBleachingPoor -0.077685 0.05687225 336 -1.36595 0.1729
## e.(Intercept) 6.917471 0.10509048 336 65.82395 0.0000
## e.ReefChilcott -0.091070 0.14820155 336 -0.61450 0.5393
## e.ReefFlinders 0.040367 0.14753203 336 0.27362 0.7845
## e.ReefHerald -0.536771 0.15546134 336 -3.45276 0.0006
## e.ReefLihou 0.057506 0.14870335 336 0.38671 0.6992
## e.ReefMoore -0.144701 0.16376371 336 -0.88359 0.3775
## Correlation:
## b d.(In) d.RfCh d.RfFl d.RfHr d.RfLh
## d.(Intercept) 0.046
## d.ReefChilcott -0.012 -0.920
## d.ReefFlinders -0.022 -0.948 0.873
## d.ReefHerald -0.018 -0.895 0.824 0.848
## d.ReefLihou -0.018 -0.930 0.856 0.882 0.832
## d.ReefMoore -0.014 -0.855 0.788 0.811 0.765 0.796
## d.catBleachingFair -0.006 -0.940 0.866 0.891 0.841 0.874
## d.catBleachingPoor -0.017 -0.946 0.871 0.897 0.847 0.880
## d.ReefChilcott:catBleachingFair 0.000 0.844 -0.910 -0.800 -0.756 -0.785
## d.ReefFlinders:catBleachingFair 0.003 0.877 -0.808 -0.920 -0.785 -0.816
## d.ReefHerald:catBleachingFair 0.003 0.816 -0.752 -0.774 -0.908 -0.760
## d.ReefLihou:catBleachingFair 0.006 0.864 -0.796 -0.820 -0.774 -0.924
## d.ReefMoore:catBleachingFair 0.000 0.795 -0.733 -0.754 -0.712 -0.740
## d.ReefChilcott:catBleachingPoor -0.006 0.738 -0.803 -0.701 -0.662 -0.688
## d.ReefFlinders:catBleachingPoor 0.015 0.758 -0.699 -0.793 -0.679 -0.706
## d.ReefHerald:catBleachingPoor 0.010 0.749 -0.690 -0.710 -0.834 -0.697
## d.ReefLihou:catBleachingPoor 0.010 0.819 -0.754 -0.776 -0.733 -0.876
## d.ReefMoore:catBleachingPoor -0.003 0.747 -0.689 -0.709 -0.670 -0.696
## e.(Intercept) -0.108 -0.126 0.112 0.117 0.110 0.114
## e.ReefChilcott 0.064 0.088 -0.154 -0.083 -0.078 -0.081
## e.ReefFlinders 0.205 0.095 -0.082 -0.150 -0.081 -0.084
## e.ReefHerald 0.047 0.084 -0.076 -0.078 -0.144 -0.077
## e.ReefLihou 0.192 0.094 -0.081 -0.085 -0.080 -0.145
## e.ReefMoore 0.131 0.083 -0.073 -0.076 -0.072 -0.074
## d.RfMr d.ctBF d.ctBP d.RC:BF d.RF:BF d.RH:BF
## d.(Intercept)
## d.ReefChilcott
## d.ReefFlinders
## d.ReefHerald
## d.ReefLihou
## d.ReefMoore
## d.catBleachingFair 0.804
## d.catBleachingPoor 0.809 0.898
## d.ReefChilcott:catBleachingFair -0.722 -0.898 -0.806
## d.ReefFlinders:catBleachingFair -0.750 -0.933 -0.838 0.838
## d.ReefHerald:catBleachingFair -0.698 -0.869 -0.780 0.780 0.811
## d.ReefLihou:catBleachingFair -0.740 -0.920 -0.826 0.826 0.859 0.799
## d.ReefMoore:catBleachingFair -0.921 -0.847 -0.760 0.761 0.790 0.736
## d.ReefChilcott:catBleachingPoor -0.632 -0.702 -0.782 0.738 0.655 0.610
## d.ReefFlinders:catBleachingPoor -0.649 -0.720 -0.802 0.647 0.742 0.625
## d.ReefHerald:catBleachingPoor -0.641 -0.711 -0.792 0.638 0.663 0.766
## d.ReefLihou:catBleachingPoor -0.700 -0.777 -0.866 0.698 0.725 0.675
## d.ReefMoore:catBleachingPoor -0.876 -0.710 -0.791 0.638 0.663 0.617
## e.(Intercept) 0.105 0.010 0.039 -0.008 -0.009 -0.008
## e.ReefChilcott -0.074 -0.007 -0.027 0.017 0.006 0.006
## e.ReefFlinders -0.076 -0.008 -0.030 0.006 0.015 0.006
## e.ReefHerald -0.070 -0.007 -0.026 0.006 0.006 0.016
## e.ReefLihou -0.076 -0.008 -0.029 0.006 0.007 0.006
## e.ReefMoore -0.166 -0.007 -0.026 0.005 0.006 0.006
## d.RL:BF d.RM:BF d.RC:BP d.RF:BP d.RH:BP d.RL:BP
## d.(Intercept)
## d.ReefChilcott
## d.ReefFlinders
## d.ReefHerald
## d.ReefLihou
## d.ReefMoore
## d.catBleachingFair
## d.catBleachingPoor
## d.ReefChilcott:catBleachingFair
## d.ReefFlinders:catBleachingFair
## d.ReefHerald:catBleachingFair
## d.ReefLihou:catBleachingFair
## d.ReefMoore:catBleachingFair 0.779
## d.ReefChilcott:catBleachingPoor 0.646 0.594
## d.ReefFlinders:catBleachingPoor 0.662 0.610 0.627
## d.ReefHerald:catBleachingPoor 0.654 0.602 0.619 0.635
## d.ReefLihou:catBleachingPoor 0.821 0.658 0.677 0.694 0.685
## d.ReefMoore:catBleachingPoor 0.653 0.814 0.619 0.634 0.626 0.685
## e.(Intercept) -0.009 -0.008 -0.028 -0.031 -0.030 -0.033
## e.ReefChilcott 0.006 0.006 0.074 0.022 0.021 0.023
## e.ReefFlinders 0.007 0.006 0.019 0.020 0.023 0.025
## e.ReefHerald 0.006 0.005 0.019 0.021 0.044 0.022
## e.ReefLihou 0.007 0.006 0.019 0.024 0.023 0.038
## e.ReefMoore 0.006 0.021 0.018 0.021 0.020 0.022
## d.RM:BP e.(In) e.RfCh e.RfFl e.RfHr e.RfLh
## d.(Intercept)
## d.ReefChilcott
## d.ReefFlinders
## d.ReefHerald
## d.ReefLihou
## d.ReefMoore
## d.catBleachingFair
## d.catBleachingPoor
## d.ReefChilcott:catBleachingFair
## d.ReefFlinders:catBleachingFair
## d.ReefHerald:catBleachingFair
## d.ReefLihou:catBleachingFair
## d.ReefMoore:catBleachingFair
## d.ReefChilcott:catBleachingPoor
## d.ReefFlinders:catBleachingPoor
## d.ReefHerald:catBleachingPoor
## d.ReefLihou:catBleachingPoor
## d.ReefMoore:catBleachingPoor
## e.(Intercept) -0.029
## e.ReefChilcott 0.020 -0.708
## e.ReefFlinders 0.020 -0.726 0.512
## e.ReefHerald 0.020 -0.673 0.477 0.486
## e.ReefLihou 0.020 -0.719 0.508 0.537 0.481
## e.ReefMoore 0.100 -0.648 0.458 0.479 0.435 0.473
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.70090770 -0.37138277 0.01415396 0.46630647 7.08204470
##
## Number of Observations: 501
## Number of Groups: 141There is a small effect of catBleaching on d but the interaction is not significant. We can see how catBleaching is affecting a subset of reefs, notably Bougainville and Moore. Bothe have small sample sizes of ‘Poor’ categories, which may explain the weak interaction. Overall, we can ignore the interaction and proceed with a simpler model that does not involve an interaction catBleaching*Reef. However, we may wish to control fo catBleaching as a random effect. Best model: reef.ahum.2
emmeans(reef.ahum.5, ~catBleaching|Reef, param = "d") %>% plot()
Furthermore, if we were to include the interaction, we have to remove some reefs where some bleaching categories are not present. We will favour including all reefs and incorporate catBleaching as a random effect.
modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = ahum.2,
curveid = Reef,
pmodels = c( ~ 1, ~ Reef, ~ Reef), bcVal = 0.5)
#random effects... can we improved the model
reef.ahum.6 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = d ~ 1|catBleaching + e ~ 1|Vial,
fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
reef.ahum.7 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = e ~ 1|catBleaching,
fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
reef.ahum.8 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = d ~ 1|catBleaching,
fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
reef.ahum.9 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = d + e ~ 1|Vial,
fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
reef.ahum.10 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = d + e ~ 1|catBleaching/Vial,
fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
reef.ahum.11 = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.2,
random = d + e ~ 1|catBleaching,
fixed = list(b ~ 1, d ~ Reef, e ~ Reef),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))
AIC(reef.ahum.2, reef.ahum.6, reef.ahum.7, reef.ahum.8, reef.ahum.9, reef.ahum.10, reef.ahum.11)## df AIC
## reef.ahum.2 15 -1221.154
## reef.ahum.6 15 -1221.154
## reef.ahum.7 15 -1221.154
## reef.ahum.8 15 -1221.154
## reef.ahum.9 17 -1217.154
## reef.ahum.10 20 -1211.154
## reef.ahum.11 17 -1217.154Based on what we saw in the species model, interaction with catBleaching above and random effects here, catBleaching may be having most of an impact on d but not e but we still want to control for Vial on e so reef.ahum.6 would be the top choice. Best: reef.ahum.6
However, the models don’t converge when consider all reefs tested with each species. So we have to simply the model further. For consistency and to compare results, we will also need to use the same model for every species.
A. cf humilis
#with all reefs
modNaive.ahum <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = ahum.1,
curveid = Reef,
pmodels = c( ~ 1, ~ 1, ~ Reef), bcVal = 0.5)
reef.ahum = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = ahum.1,
random = d ~ 1|catBleaching + e ~ 1|Vial,
fixed = list(b ~ 1, d ~ 1, e ~ Reef),
start = coef(modNaive.ahum),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))To protect our original models, we’ve saved them as a Rdata file
load("3_outputs/nlme_A.humilis by reef.RData")
emmeans(reef.ahum, ~Reef, param = "e") %>% plot(comparisons = TRUE)
emmeans(reef.ahum, pairwise~Reef, param = "e")$contrast %>% plot()
emmeans(reef.ahum, ~Reef, param = 'e', level = 0.95) ## Reef emmean SE df lower.CL upper.CL
## Bougainville 6.70 0.0941 468 6.52 6.89
## Chilcott 6.91 0.1107 468 6.70 7.13
## Flinders 6.93 0.0986 468 6.74 7.13
## Frederick 6.80 0.1194 468 6.56 7.03
## Herald 6.37 0.1156 468 6.14 6.59
## Lihou 7.01 0.1050 468 6.80 7.22
## Moore 6.64 0.1112 468 6.42 6.86
## Wreck 8.26 0.0885 468 8.08 8.43
##
## Confidence level used: 0.95emmeans(reef.ahum, pairwise~Reef, param = "e")$contrast## contrast estimate SE df t.ratio p.value
## Bougainville - Chilcott -0.2082 0.144 468 -1.447 0.8348
## Bougainville - Flinders -0.2268 0.135 468 -1.677 0.7024
## Bougainville - Frederick -0.0941 0.151 468 -0.625 0.9985
## Bougainville - Herald 0.3386 0.148 468 2.288 0.3021
## Bougainville - Lihou -0.3054 0.140 468 -2.186 0.3624
## Bougainville - Moore 0.0617 0.145 468 0.426 0.9999
## Bougainville - Wreck -1.5534 0.127 468 -12.211 <.0001
## Chilcott - Flinders -0.0186 0.147 468 -0.127 1.0000
## Chilcott - Frederick 0.1141 0.161 468 0.708 0.9968
## Chilcott - Herald 0.5468 0.159 468 3.441 0.0145
## Chilcott - Lihou -0.0972 0.151 468 -0.644 0.9982
## Chilcott - Moore 0.2699 0.156 468 1.733 0.6654
## Chilcott - Wreck -1.3452 0.140 468 -9.594 <.0001
## Flinders - Frederick 0.1327 0.153 468 0.869 0.9886
## Flinders - Herald 0.5655 0.151 468 3.746 0.0049
## Flinders - Lihou -0.0786 0.139 468 -0.565 0.9992
## Flinders - Moore 0.2885 0.145 468 1.984 0.4942
## Flinders - Wreck -1.3266 0.135 468 -9.853 <.0001
## Frederick - Herald 0.4328 0.165 468 2.624 0.1499
## Frederick - Lihou -0.2113 0.157 468 -1.349 0.8794
## Frederick - Moore 0.1558 0.162 468 0.965 0.9791
## Frederick - Wreck -1.4592 0.147 468 -9.899 <.0001
## Herald - Lihou -0.6440 0.155 468 -4.155 0.0010
## Herald - Moore -0.2769 0.160 468 -1.735 0.6641
## Herald - Wreck -1.8920 0.144 468 -13.117 <.0001
## Lihou - Moore 0.3671 0.150 468 2.451 0.2194
## Lihou - Wreck -1.2480 0.139 468 -8.982 <.0001
## Moore - Wreck -1.6151 0.143 468 -11.299 <.0001
##
## P value adjustment: tukey method for comparing a family of 8 estimates#Some significant differencesP. meandrina
#with all reefs
modNaive.pmea <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = pmea.1,
curveid = Reef,
pmodels = c( ~ 1, ~ 1, ~ Reef), bcVal = 0.5)
reef.pmea = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = pmea.1,
random = d ~ 1|catBleaching + e ~ 1|Vial,
fixed = list(b ~ 1, d ~ 1, e ~ Reef),
start = coef(modNaive.pmea),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))load("3_outputs/nlme_P.meandrina by reef.RData")
emmeans(reef.pmea, ~Reef, param = "e") %>% plot(comparisons = TRUE)
emmeans(reef.pmea, pairwise~Reef, param = "e")$contrast %>% plot()
emmeans(reef.pmea, ~Reef, param = 'e', level = 0.95) ## Reef emmean SE df lower.CL upper.CL
## Bougainville 7.09 0.0973 271 6.89 7.28
## Chilcott 7.08 0.1148 271 6.85 7.31
## Flinders 8.11 0.1555 271 7.81 8.42
## Frederick 7.06 0.1166 271 6.83 7.29
## Herald 7.25 0.1885 271 6.88 7.62
## Lihou 6.96 0.1770 271 6.61 7.31
## Moore 7.26 0.1159 271 7.03 7.49
## Saumarez 7.86 0.1167 271 7.63 8.09
## Wreck 8.03 0.1016 271 7.83 8.23
##
## Confidence level used: 0.95emmeans(reef.pmea, pairwise~Reef, param = "e")$contrast## contrast estimate SE df t.ratio p.value
## Bougainville - Chilcott 0.00496 0.150 271 0.033 1.0000
## Bougainville - Flinders -1.02701 0.180 271 -5.717 <.0001
## Bougainville - Frederick 0.02894 0.148 271 0.195 1.0000
## Bougainville - Herald -0.16248 0.211 271 -0.771 0.9975
## Bougainville - Lihou 0.12858 0.203 271 0.633 0.9994
## Bougainville - Moore -0.17636 0.149 271 -1.184 0.9593
## Bougainville - Saumarez -0.77570 0.146 271 -5.298 <.0001
## Bougainville - Wreck -0.94347 0.136 271 -6.951 <.0001
## Chilcott - Flinders -1.03197 0.194 271 -5.318 <.0001
## Chilcott - Frederick 0.02398 0.163 271 0.147 1.0000
## Chilcott - Herald -0.16743 0.220 271 -0.762 0.9977
## Chilcott - Lihou 0.12362 0.208 271 0.594 0.9996
## Chilcott - Moore -0.18131 0.161 271 -1.123 0.9704
## Chilcott - Saumarez -0.78065 0.164 271 -4.758 0.0001
## Chilcott - Wreck -0.94842 0.154 271 -6.144 <.0001
## Flinders - Frederick 1.05595 0.188 271 5.603 <.0001
## Flinders - Herald 0.86454 0.243 271 3.562 0.0128
## Flinders - Lihou 1.15559 0.241 271 4.792 0.0001
## Flinders - Moore 0.85065 0.191 271 4.445 0.0004
## Flinders - Saumarez 0.25131 0.184 271 1.366 0.9095
## Flinders - Wreck 0.08355 0.176 271 0.475 0.9999
## Frederick - Herald -0.19141 0.220 271 -0.871 0.9943
## Frederick - Lihou 0.09964 0.214 271 0.465 0.9999
## Frederick - Moore -0.20529 0.161 271 -1.273 0.9382
## Frederick - Saumarez -0.80463 0.157 271 -5.132 <.0001
## Frederick - Wreck -0.97240 0.147 271 -6.605 <.0001
## Herald - Lihou 0.29105 0.258 271 1.128 0.9695
## Herald - Moore -0.01388 0.219 271 -0.063 1.0000
## Herald - Saumarez -0.61322 0.219 271 -2.797 0.1212
## Herald - Wreck -0.78099 0.212 271 -3.678 0.0085
## Lihou - Moore -0.30493 0.211 271 -1.446 0.8788
## Lihou - Saumarez -0.90427 0.218 271 -4.147 0.0015
## Lihou - Wreck -1.07204 0.211 271 -5.090 <.0001
## Moore - Saumarez -0.59934 0.160 271 -3.737 0.0069
## Moore - Wreck -0.76711 0.151 271 -5.082 <.0001
## Saumarez - Wreck -0.16777 0.141 271 -1.187 0.9587
##
## P value adjustment: tukey method for comparing a family of 9 estimates#Some significant differencesP. verrucosa
#with all reefs
modNaive.pver <- drm(median.yield ~ MMM_dev, fct = L.3(),
data = pver.1,
curveid = Reef,
pmodels = c( ~ 1, ~ 1, ~ Reef), bcVal = 0.5)
reef.pver = nlme(median.yield ~ NLS.L3(MMM_dev, b, d, e), data = pver.1,
random = d ~ 1|catBleaching + e ~ 1|Vial,
fixed = list(b ~ 1, d ~ 1, e ~ Reef),
start = coef(modNaive.pver),
control = list(msMaxIter = 2000, maxIter=2000, minScale=0.0001))load("3_outputs/nlme_P.verrucosa by reef.RData")
emmeans(reef.pver, ~Reef, param = "e") %>% plot(comparisons = TRUE)
emmeans(reef.pver, pairwise~Reef, param = "e")$contrast %>% plot()
emmeans(reef.pver, ~Reef, param = 'e', level = 0.95) ## Reef emmean SE df lower.CL upper.CL
## Bougainville 7.51 0.0587 229 7.40 7.63
## Flinders 8.11 0.0903 229 7.93 8.28
## Herald 7.26 0.1251 229 7.01 7.51
## Lihou 7.74 0.1307 229 7.48 7.99
## Moore 7.54 0.1373 229 7.27 7.81
## Wreck 8.09 0.1279 229 7.83 8.34
##
## Confidence level used: 0.95emmeans(reef.pver, pairwise~Reef, param = "e")$contrast## contrast estimate SE df t.ratio p.value
## Bougainville - Flinders -0.5949 0.0938 229 -6.339 <.0001
## Bougainville - Herald 0.2514 0.1373 229 1.830 0.4483
## Bougainville - Lihou -0.2240 0.1311 229 -1.709 0.5275
## Bougainville - Moore -0.0317 0.1418 229 -0.224 0.9999
## Bougainville - Wreck -0.5736 0.1311 229 -4.375 0.0003
## Flinders - Herald 0.8463 0.1536 229 5.511 <.0001
## Flinders - Lihou 0.3708 0.1314 229 2.822 0.0576
## Flinders - Moore 0.5631 0.1481 229 3.801 0.0025
## Flinders - Wreck 0.0212 0.1348 229 0.157 1.0000
## Herald - Lihou -0.4754 0.1798 229 -2.644 0.0910
## Herald - Moore -0.2832 0.1844 229 -1.535 0.6418
## Herald - Wreck -0.8250 0.1785 229 -4.623 0.0001
## Lihou - Moore 0.1923 0.1728 229 1.112 0.8759
## Lihou - Wreck -0.3496 0.1611 229 -2.170 0.2559
## Moore - Wreck -0.5419 0.1746 229 -3.104 0.0259
##
## P value adjustment: tukey method for comparing a family of 6 estimates#Some significant differencesAbsolute Temperature ED50s (Table S5 SOM)
ahum.abs = ahum.1 %>% group_by(Reef) %>%
summarise(MMM = mean(MMM)) %>%
left_join(emmeans(reef.ahum, ~Reef, param = "e") %>% as.data.frame(), by = "Reef") %>%
mutate(absolute = MMM + emmean) %>%
mutate(absolute.lowerCL = MMM + lower.CL) %>%
mutate(absolute.upperCL = MMM + upper.CL)
pmea.abs = pmea.1%>% group_by(Reef) %>%
summarise(MMM = mean(MMM)) %>%
left_join(emmeans(reef.pmea, ~Reef, param = "e") %>% as.data.frame(), by = "Reef") %>%
mutate(absolute = MMM + emmean) %>%
mutate(absolute.lowerCL = MMM + lower.CL) %>%
mutate(absolute.upperCL = MMM + upper.CL)
pver.abs = pver.1 %>% group_by(Reef) %>%
summarise(MMM = mean(MMM)) %>%
left_join(emmeans(reef.pver, ~Reef, param = "e") %>% as.data.frame(), by = "Reef") %>%
mutate(absolute = MMM + emmean) %>%
mutate(absolute.lowerCL = MMM + lower.CL) %>%
mutate(absolute.upperCL = MMM + upper.CL)2.5 Model predictions: Species*Reef
##Create predictions data frame
predframe.ahum = with(ahum.1, expand.grid(Reef = levels(Reef), MMM_dev=seq(min(MMM_dev), max(MMM_dev) + 1, 0.1)))
#Adding confidence intervals
prdc.ahum = as.data.frame(nlraa::predict_nlme(reef.ahum, newdata = predframe.ahum, interval = "confidence", level = 0.95))
prdc.ahum = bind_cols(predframe.ahum, prdc.ahum)
# Predict ED50 for each species
#(Since we're not varyin d, we don't have to do this but could be handy nevertheless. Estimate values should not change much between reefs)
ED50.ahum = emmeans(reef.ahum, ~Reef, param= "e") %>% as_tibble() %>%
dplyr::select(Reef, emmean) %>%
dplyr::rename(MMM_dev = emmean)
prdc.ED50.ahum = as.data.frame(nlraa::predict_nlme(reef.ahum, newdata = ED50.ahum, interval = "confidence", level = 0.95))
prdc.ED50.ahum = bind_cols(ED50.ahum, prdc.ED50.ahum)P. meandrina
##Create predictions data frame
predframe.pmea = with(pmea.1, expand.grid(Reef = levels(Reef), MMM_dev=seq(min(MMM_dev), max(MMM_dev) + 1, 0.1)))
#Adding confidence intervals
prdc.pmea = as.data.frame(nlraa::predict_nlme(reef.pmea, newdata = predframe.pmea, interval = "confidence", level = 0.95))
prdc.pmea = bind_cols(predframe.pmea, prdc.pmea)
# Predict ED50 for each species
ED50.pmea = emmeans(reef.pmea, ~Reef, param= "e") %>% as_tibble() %>%
dplyr::select(Reef, emmean) %>%
dplyr::rename(MMM_dev = emmean)
prdc.ED50.pmea = as.data.frame(nlraa::predict_nlme(reef.pmea, newdata = ED50.pmea, interval = "confidence", level = 0.95))
prdc.ED50.pmea = bind_cols(ED50.pmea, prdc.ED50.pmea)P. verrucosa
##Create predictions data frame
predframe.pver = with(pver.1,
expand.grid(Reef = levels(Reef),
MMM_dev=seq(min(MMM_dev), max(MMM_dev) + 1, 0.1)))
#Adding confidence intervals
prdc.pver = as.data.frame(nlraa::predict_nlme(reef.pver, newdata = predframe.pver,
interval = "confidence", level = 0.95))
prdc.pver = bind_cols(predframe.pver, prdc.pver)
# Predict ED50 for each species
ED50.pver = emmeans(reef.pver, ~Reef, param= "e") %>% as_tibble() %>%
dplyr::select(Reef, emmean) %>%
dplyr::rename(MMM_dev = emmean)
prdc.ED50.pver = as.data.frame(nlraa::predict_nlme(reef.pver, newdata = ED50.pver, interval = "confidence", level = 0.95))
prdc.ED50.pver = bind_cols(ED50.pver, prdc.ED50.pver)2.6 Figure: Reef comparison
#Base order of reefs and colours on the average ED50.
Reef.ED50 = prdc.ED50.ahum %>% mutate(Species = "A. cf humilis") %>% dplyr::select(Species, Reef, MMM_dev) %>%
bind_rows(prdc.ED50.pmea %>% mutate(Species = "P. meandrina") %>% dplyr::select(Species, Reef, MMM_dev)) %>%
bind_rows(prdc.ED50.pver %>% mutate(Species = "P. verrucosa") %>% dplyr::select(Species, Reef, MMM_dev)) %>%
spread(key = Species, value = MMM_dev) %>%
rowwise(Reef) %>%
dplyr::summarise(average = mean(c(`A. cf humilis`, `P. meandrina`, `P. verrucosa`), na.rm = T)) %>%
arrange(average) ## `summarise()` has grouped output by 'Reef'. You can override using the
## `.groups` argument.Reef.order = fct_reorder(Reef.ED50$Reef, Reef.ED50$average)
#Combine all data and reorder Reef factor
line.pred.dat = prdc.ahum %>% mutate(Species = "A. cf humilis") %>%
bind_rows(prdc.pmea %>% mutate(Species = "P. meandrina")) %>%
bind_rows(prdc.pver %>% mutate(Species = "P. verrucosa")) %>%
mutate(Reef = factor(Reef, levels(Reef.order)))
point.obs.dat = ahum.1 %>%
bind_rows(pmea.1) %>%
bind_rows(pver.1) %>%
mutate(Reef = factor(Reef, levels(Reef.order)))
point.ED50.dat = prdc.ED50.ahum %>% mutate(Species = "A. cf humilis") %>%
bind_rows(prdc.ED50.pmea %>% mutate(Species = "P. meandrina")) %>%
bind_rows(prdc.ED50.pver %>% mutate(Species = "P. verrucosa")) %>%
mutate(Reef = factor(Reef, levels(Reef.order)))Fig 2B: PSII comparison by Reef
col = wes_palette("Zissou1", 9, type = "continuous")
Reef.plot = ggplot() +
geom_segment(data = point.ED50.dat, aes(x = MMM_dev, xend = MMM_dev, y = -Inf, yend = Estimate , color = Reef),
linetype = 1, size = .3) +
geom_line(data = line.pred.dat, aes(colour=Reef, y=Estimate, x=MMM_dev), size = .3) +
geom_point(data = point.ED50.dat, aes(x = MMM_dev, y = Estimate, fill = Reef), size = 1.5, stroke = .2, shape = 21, col = "white") +
geom_text(data = data.frame(Species = c("A. cf humilis", "P. meandrina", "P. verrucosa")),
aes(y = .71, x = 2, label = Species), size = 3, fontface = "italic", family = "Helvetica", hjust = 0) +
#scale_y_continuous(labels = scales::percent, breaks = c(0, .5, 1)) +
scale_x_continuous(breaks = c(3,6,9)) +
coord_cartesian(xlim = c(2,9.6), ylim = c(0,.75)) +
scale_fill_manual(values=col) +
scale_colour_manual(values=col) +
facet_grid(Species~.) +
labs(x="Temp. above local MMM (°C) ",
y="PSII-yield") +
theme_classic() +
theme(axis.line = element_blank(),
panel.border = element_rect(size = .5, fill = "transparent"),
strip.text = element_blank(),
strip.background = element_blank(),
legend.title = element_blank(),
legend.background = element_rect(fill = "transparent"),
legend.position = "right",
legend.text = element_text(size = 7, family = "Helvetica"),
legend.key.size = unit(2, units = "mm"),
legend.margin = margin(0,0,0,0, unit = "mm"),
axis.text.y = element_text(size = 8, family = "Helvetica"),
axis.title = element_text(size = 9, family = "Helvetica"))
Reef.plot
2.6 Model comparison
Does it make a difference to estimate e using the 3 way interaction of for each reef separately?
Fig SX: model comparison
#All species in the same model
mod.all = emmeans(reef.1, ~Reef|Species, param = "e") %>% as.tibble() %>%
dplyr::select(Species, Reef, emmean, lower.CL, upper.CL)## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.#Species modelled separately
mod.sep = emmeans(reef.ahum, ~Reef, param= "e") %>% as.tibble() %>%
mutate(Species = "A. cf humilis") %>% dplyr::select(Species, Reef, emmean, lower.CL, upper.CL) %>%
bind_rows(emmeans(reef.pmea, ~Reef, param= "e") %>% as.tibble() %>%
mutate(Species = "P. meandrina") %>% dplyr::select(Species, Reef, emmean, lower.CL, upper.CL)) %>%
bind_rows(emmeans(reef.pver, ~Reef, param= "e") %>% as.tibble() %>%
mutate(Species = "P. verrucosa") %>% dplyr::select(Species, Reef, emmean, lower.CL, upper.CL)) %>%
mutate(Reef = factor(Reef, levels(Reef.order)))
ggplot(bind_rows(mod.all %>% mutate(model = "combined"), mod.sep %>% mutate(model = "separate")),
aes(y = emmean, x = Reef, ymin = lower.CL, ymax = upper.CL, col = model)) +
geom_pointrange(position = position_dodge(width = .5)) +
scale_color_manual(values = c("grey30", "skyblue2")) +
facet_wrap(~Species, scales = "free_x") +
labs(y = "ED50 estimate", x = "") +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.line = element_blank(),
panel.border = element_rect(color = "black", fill = "transparent", size = 1),
strip.background = element_blank(),
strip.text = element_text(face = "italic"),
legend.title = element_blank(),
legend.position = c(.9,.15),
legend.background = element_rect(fill = "transparent"))
ggsave("2_figures/FigS4_Model comparison.png", width = 140, height = 70, units = "mm")It makes no difference what model we use.
3. DRC Results
Table S1. ED50 values for each reef nested within species + 95 CI
library(gt)
library(tidyverse)
tableED50reef = mod.sep %>%
mutate(emmean = round(emmean,2), CI = paste(round(lower.CL,2), round(upper.CL,2), sep = "-")) %>%
dplyr::select(-lower.CL, -upper.CL) %>%
mutate(ED50 = paste(emmean, CI, sep = " ")) %>%
dplyr::select(-emmean, -CI) %>%
spread(key = Species, value = ED50) %>%
mutate(Reef = factor(Reef, levels = Reef.order)) %>% arrange(Reef)
gt_tbl.reefED50 <- gt(tableED50reef) %>%
tab_header(
title = "Table S4. Summary tables derived from emmeans for ED50 values related to each reef and species combination (95% confidence intervals). ",) %>%
tab_spanner(label = "Species",
columns = c("A. cf humilis", "P. meandrina", "P. verrucosa"))
#gt_tbl.reefED50Figure 2
g1 = ggarrange(
#Species plot
Species.plot + labs(x = "") +
theme(axis.text.x = element_blank(),
plot.margin = margin(0.1,0.1,0,.6, unit = "cm")),
#Reef plot
Reef.plot + theme(plot.margin = margin(0.1,0.24,0.1,.6, unit = "cm")),
nrow = 2, heights = c(.33,.66), labels = c("a", "b"),
font.label = list(size = 12, family = "Helvetica"))
g1
ggsave("2_figures/Fig2_v1.png", dpi =300, width=90, height=110, units = "mm")#4. Drivers of heat tolerance among reefs in the CSMP
There are significant differences in the estimated ED50 among reefs. Can different climatic regimes be driving these patterns?
4.2 Data
Fig 3a . PSII-50 by reef and species with 95%CI
col = wes_palette("GrandBudapest1", 4, type = "continuous")
Fig3a = ggplot(mod.sep, aes(y = emmean, x = Reef, ymin = lower.CL, ymax = upper.CL, col = Species)) +
geom_pointrange(position = position_dodge(width = .4), size = .2) +
scale_y_continuous(expand = c(0,0), lim = c(5.9,8.5), breaks = c(6.5, 7, 7.5, 8)) +
scale_colour_manual(values=col) +
labs(x="",
y="ED50 (95% CI)") +
theme_classic() +
theme(axis.line = element_blank(),
panel.border = element_rect(size = .5, fill = "transparent"),
legend.title = element_blank(),
legend.background = element_rect(fill = "transparent"),
legend.position = c(.22,.9),
legend.text = element_text(face = "italic", size = 10, family = "Helvetica"),
legend.key.size = unit(3, units = "mm"),
legend.margin = margin(0,0,0,0, unit = "mm"),
axis.text = element_text(size = 14, family = "Helvetica"),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.title = element_text(size = 14, family = "Helvetica"))
Fig3a
We can observe a large variation in the point of inflection between
reefs in the CSMP that span several degrees of latitude and an array of
climatic regimes.
Climate predictors of PSII-50
We will use the PSII-50 values derived from separate reef models and test against a suite of environmental variables
Thermal history metadata (based on NOAA ….. ) DHW3: Number of Degree Heating Weeks (DHW) events greater than 3 between 1985 - 2020 DHW4: Number of DHW events greater than 4C-weeks between 1985 - 2020 DHW6: Number of DHW events greater than 6C-weeks between 1985 - 2020 DHW8: Number of DHW events greater than 8C-weeks between 1985 - 2020 DHW9: Number of DHW events greater than 9C-weeks between 1985 - 2020 maxDHW: Maximum DHW experienced between 1985 - 2020 meanDHW: Average DHWs experienced between 1985 - 2020
MMM: Maximum Monthly Mean SST between 1986 - 2010 meanSST: Average Sea Surface Temperature (SST) between 1985 - 2020 maxSST: Maximum Sea Surface Temperature (SST) between 1985 - 2020 minSST: Max Sea Surface Temperature (SST) between 1985 - 2020 rangeSST: Range between Max and Min SST values between 1985 - 2020 varSST: Variability in SST between 1985 - 2020
Recent.maxDHW Maximum DHW exposue between 2016 - 2020 Recent.meanDHW Average maximum DHW exposure between 2016 - 2020
ReturnDHW3 The return time between events where DHW exceeded 3C-weeks ReturnDHW4 The return time between events where DHW exceeded 4C-weeks
Merging all data frames together
#Load data for thermal history data
load("1_data/SiteDisturbanceHistory_DHW.RData")
load("1_data/reef.info.RData")
#DHW values for all reefs at the time of sampling (Feb-Mar 2020)
DHW2020 = bind_rows(ahum.1, pmea.1, pver.1) %>%
group_by(Reef, Site, DHW) %>%
dplyr::summarise() %>%
group_by(Reef) %>%
dplyr::summarise(DHW2020 = mean(DHW, na.rm = TRUE)) ## `summarise()` has grouped output by 'Reef', 'Site'. You can override using the
## `.groups` argument.#we merge data frames for thermal history metrics and geographic information with the ED50 values
clim.variables = site.bleachings %>% ungroup() %>% dplyr::select(-Sector, -Site, -Month) %>%
group_by(Reef) %>% dplyr::summarise_all(mean) %>%
left_join(DHW2020, by = "Reef") %>%
left_join(reef.info %>% ungroup() %>% dplyr::select(-SectorRegion, -Sector, -Region), by = "Reef") %>%
filter(Reef %in% mod.sep$Reef)
psii.prediction = clim.variables %>%
gather(key = variable, value = value, -Reef) %>%
full_join(mod.sep %>% dplyr::rename(ED50 = emmean), by = "Reef") %>%
filter(! is.na(value)) %>%
dplyr::mutate(value = round(value, 2))
#linear relationships
col = wes_palette("GrandBudapest1", 4, type = "continuous")
ggplot(psii.prediction, aes(y = ED50, x = value, col = Species, fill = Species)) +
geom_point() +
geom_smooth(method = "lm") +
scale_fill_manual(values=col) +
scale_colour_manual(values=col) +
labs(y = "ED50", x = "Environmental metric") +
facet_wrap(~variable, scales = "free") +
theme_classic() +
theme(axis.line = element_blank(),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.text = element_text(size = 8, family = "Helvetica"),
strip.text = element_text(size = 7)) ## `geom_smooth()` using formula 'y ~ x'
ggsave("2_figures/FigS5_correlation plot.png", height = 5, width = 7, units = "in")## `geom_smooth()` using formula 'y ~ x'At first glance, we see that many variables show some relationship with ED50, but that heat tolerance drivers might be species-specific (ED50) in response to environmental variables. For example, DHW3 has strong correlation to ED50 value for A. humilis but the trend looks less strong for either of the Poci species, which may point towards an interaction.
4.1 Correlation matrix of thermal history metrics
Correlations between climate variables
Before we interpret these relationships, it’s important to consider that many of these variables are co-linear. Any metrics that have high correlation may represent the same climate variability and do not need to be included.
library(corrplot)## corrplot 0.92 loadedclim.variables %>% dplyr::select(-Reef) %>% cor() %>% corrplot()
Correlations between ED50 values and thermal variables
Let’s look at correlation coefficients from these data to see which variables show strong correlation.
psii.prediction %>%
group_by(variable, Species) %>%
summarise(correlation = cor(value, ED50)) %>%
ggplot(aes(y = variable, x = Species, fill = correlation)) +
geom_tile() +
scale_fill_gradient2(low = "darkred", mid = "white", high = "blue3") +
theme_classic() + theme(axis.line = element_blank())
ggsave("2_figures/FigS7_correlation by species.png", height = 5, width = 7, units = "in")
save(psii.prediction, file = "1_data/psii.prediction.Rdata")Looking at each species correlation values, we find that each species responds differently to each environmental variables. In common, all species seem to have a very fairly strong positive relationship to varSST, and negative relationship to recent meanDHW.
Here we are looking at the correlation values between enviro and ED50 selecting a subset with any variables that are > 0.40.
psii.prediction %>%
group_by(variable) %>%
summarise(correlation = cor(value, ED50))## # A tibble: 24 × 2
## variable correlation
## <chr> <dbl>
## 1 Complexity -0.202
## 2 DHW2 0.353
## 3 DHW2020 -0.501
## 4 DHW3 0.353
## 5 DHW4 0.458
## 6 DHW6 -0.0234
## 7 DHW8 -0.375
## 8 DHW9 -0.177
## 9 Lat -0.498
## 10 Long 0.365
## # … with 14 more rows
## # ℹ Use `print(n = ...)` to see more rowspsii.wide = psii.prediction %>% spread(variable, value)
save(psii.wide, file = "1_data/psii.wide.Rdata")We can take a subset of variables to represent the range of climate regimes and that have low colinearity to other variables. We are excluding MMM and DHW2020 because the collinearity is > 0.80 to other variables that are of importance and have a stronger correlation to ED50s (e.g. recent.meanDHW).
climate.subset = c("DHW2020", "DHW4", "Lat", "meanSST", "minSST", "MMM", "rangeSST", "recent.meanDHW", "returnDHW6", "varSST")
clim.variables %>% dplyr::select(all_of(climate.subset)) %>%
cor() %>% corrplot() 
#we will keep all but range SST. We are now subsetting variables for any that have high correlation to ED50. - First, we excluded latitude, meanSST, minSST, maxSST because they are all highly correlated to each other. - Then, we also exclude DHW2020 and MMM because they are highly correlated to recent.meanDHW. We choose to keep recent.meanDHW because it has the highest correlation with ED50 out of any of these 3 variables. - We then exclude rangeSST because of high correlation with varSST. VarSST has stronger relationship with ED50 so we preferentially keep it in. This leaves us with 4 remaining variables to go into the first dredge model: DHW4, recent.meanDHW, returnDHW6, and varSST.
thermal.subset = c("DHW4", "recent.meanDHW", "returnDHW6", "varSST")
clim.variables%>% dplyr::select(all_of(thermal.subset)) %>%
cor() %>% corrplot()
Below, the four plots represent some of the strongest climate variables that are not too strongly correlated with each other, but show some correlation with ED50. This allows us to assess whether some variables might have an interaction with species. We can use this subset to investigate the drivers of ED50 in a dredge model.
col = wes_palette("GrandBudapest1", 4, type = "continuous")
psii.prediction %>% filter(variable %in% thermal.subset) %>%
ggplot(aes(y = ED50, x = value, col = Species, fill = Species)) +
geom_smooth(method = "lm") +
geom_point() +
scale_fill_manual(values=col) +
scale_colour_manual(values=col) +
labs(y = "ED50") +
facet_wrap(~variable, scales = "free") +
theme_classic() +
theme(axis.line = element_blank(),
panel.border = element_rect(size = .5, fill = "transparent"),
axis.text = element_text(size = 8, family = "Helvetica"),
axis.title = element_text(size = 9, family = "Helvetica"))## `geom_smooth()` using formula 'y ~ x'
ggsave("2_figures/FigS8_EnviroRelationships.png", height = 5, width = 7, units = "in")## `geom_smooth()` using formula 'y ~ x'4.3 Dredging: Candidate Model Selection
Initial dredge
library(MuMIn)
library(lme4)## Warning: package 'lme4' was built under R version 4.0.5## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmListoptions(na.action = "na.fail")
#creating a data frame with only the 4 variables selected (DHW4, recent.meanDHW, returnDHW6, and varSST)
dredge.dat1 = psii.prediction %>% filter(variable %in% thermal.subset) %>%
spread(key = variable, value = value) %>%
dplyr::select(-ends_with(".CL"), -Reef)
#dredge model that accounts for species + the four variables in the model.
dredge.mod1 <- dredge(lm(ED50~Species +., data = dredge.dat1), rank = "AIC", m.lim = c(1,4))## Fixed term is "(Intercept)"#just look at the top 10 scoring models
head(dredge.mod1, 10)## Global model call: lm(formula = ED50 ~ Species + ., data = dredge.dat1)
## ---
## Model selection table
## (Int) DHW rcn.mDH rDH Spc vSS df logLik AIC delta weight
## 16 6.607 0.12220 -0.1630 0.05856 + 7 -0.548 15.1 0.00 0.471
## 14 5.212 0.17540 0.08732 + 6 -2.592 17.2 2.09 0.166
## 30 4.335 0.14940 0.06693 + 0.4729 7 -2.277 18.6 3.46 0.084
## 12 7.998 0.06906 -0.2910 + 6 -3.664 19.3 4.23 0.057
## 28 5.601 0.06648 -0.1862 + 0.7482 7 -2.720 19.4 4.34 0.054
## 11 8.734 -0.3401 + 5 -4.907 19.8 4.72 0.044
## 27 6.183 -0.2278 + 0.7871 6 -3.967 19.9 4.84 0.042
## 15 8.330 -0.2959 0.02869 + 6 -4.125 20.2 5.15 0.036
## 26 2.637 0.08607 + 1.4910 6 -4.335 20.7 5.57 0.029
## 31 6.697 -0.2356 0.01616 + 0.5583 7 -3.789 21.6 6.48 0.018
## Models ranked by AIC(x)#extracting the best fit model using AIC score
bestmodel <- get.models(dredge.mod1, 1)[[1]]
lm.dredge1 <- lm(bestmodel, data = dredge.dat1)
summary(lm.dredge1)##
## Call:
## lm(formula = bestmodel, data = dredge.dat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.44268 -0.15135 0.02352 0.12432 0.53704
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.60656 0.83588 7.904 4.3e-07 ***
## DHW4 0.12225 0.04909 2.490 0.023399 *
## recent.meanDHW -0.16296 0.08962 -1.818 0.086680 .
## returnDHW6 0.05856 0.02546 2.300 0.034371 *
## SpeciesP. meandrina 0.41367 0.14079 2.938 0.009184 **
## SpeciesP. verrucosa 0.68668 0.16006 4.290 0.000495 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2882 on 17 degrees of freedom
## Multiple R-squared: 0.7835, Adjusted R-squared: 0.7199
## F-statistic: 12.31 on 5 and 17 DF, p-value: 3.718e-05ggeffects::ggeffect(lm.dredge1) %>% plot## $DHW4
##
## $recent.meanDHW
##
## $returnDHW6
##
## $Species
The best model is ED50 ~ DHW4 + recent.meanDHW + returnDHW6 +
Species Only varSST was not included. However, we need to also
consider if the species interaction is significant.
Dredging with species interaction.
Dredging across all selected variables suggest a model that includes DHW4, MMM and recent.meanDHW best represents the data so far. However, we’ve seen there may be an interaction between Species * DHW4. Before forcing this interaction using dredge, we will see which model yields best AIC score with interaction and additional fixed effects.
lm.1 = lm(ED50~Species + DHW4 + recent.meanDHW + returnDHW6, data = dredge.dat1)
lm.2 = lm(ED50~Species * DHW4 + recent.meanDHW + returnDHW6, data = dredge.dat1)
lm.3 = lm(ED50~Species * returnDHW6 + DHW4 + recent.meanDHW, data = dredge.dat1)
lm.4 = lm(ED50~Species * DHW4 * returnDHW6 + recent.meanDHW, data = dredge.dat1)
AICc(lm.1, lm.2, lm.3, lm.4)## df AICc
## lm.1 7 22.56228
## lm.2 9 21.94864
## lm.3 9 28.24851
## lm.4 14 61.85056#lm.2 has the lowest AIC, suggesting it is the better model.
anova(lm.1, lm.2)## Analysis of Variance Table
##
## Model 1: ED50 ~ Species + DHW4 + recent.meanDHW + returnDHW6
## Model 2: ED50 ~ Species * DHW4 + recent.meanDHW + returnDHW6
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 17 1.41235
## 2 15 0.87572 2 0.53663 4.5959 0.02774 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#The two models are also sufficiently different from each other to include the interaction
summary(lm.2)##
## Call:
## lm(formula = ED50 ~ Species * DHW4 + recent.meanDHW + returnDHW6,
## data = dredge.dat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38367 -0.08972 0.03304 0.09462 0.43907
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.78952 0.76071 7.611 1.58e-06 ***
## SpeciesP. meandrina 1.84012 0.55064 3.342 0.004460 **
## SpeciesP. verrucosa 2.19242 0.60783 3.607 0.002589 **
## DHW4 0.25490 0.06028 4.229 0.000729 ***
## recent.meanDHW -0.17568 0.07644 -2.298 0.036344 *
## returnDHW6 0.05255 0.02162 2.431 0.028076 *
## SpeciesP. meandrina:DHW4 -0.20679 0.07812 -2.647 0.018295 *
## SpeciesP. verrucosa:DHW4 -0.21270 0.08280 -2.569 0.021379 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2416 on 15 degrees of freedom
## Multiple R-squared: 0.8658, Adjusted R-squared: 0.8031
## F-statistic: 13.82 on 7 and 15 DF, p-value: 1.633e-05summary(lm.3)##
## Call:
## lm(formula = ED50 ~ Species * returnDHW6 + DHW4 + recent.meanDHW,
## data = dredge.dat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.43125 -0.12992 0.00647 0.14245 0.57188
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.71358 0.80810 8.308 5.39e-07 ***
## SpeciesP. meandrina -0.08909 0.30582 -0.291 0.7748
## SpeciesP. verrucosa 0.35133 0.37683 0.932 0.3659
## returnDHW6 0.01414 0.03497 0.404 0.6917
## DHW4 0.13149 0.04793 2.744 0.0151 *
## recent.meanDHW -0.14462 0.08738 -1.655 0.1187
## SpeciesP. meandrina:returnDHW6 0.08080 0.04397 1.838 0.0860 .
## SpeciesP. verrucosa:returnDHW6 0.05375 0.05215 1.031 0.3190
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2771 on 15 degrees of freedom
## Multiple R-squared: 0.8235, Adjusted R-squared: 0.7411
## F-statistic: 9.997 on 7 and 15 DF, p-value: 0.0001144#the interaction with DHW4 is important but not the interaction with returnDHW6 so will be used in the full dredge model.The interaction is significant and the two models are significantly different from one another so we should include the interaction. (We also explored a three-way interaction with MMM and the interaction between Species and MMM, neither were significant)
vif(lm.2)## GVIF Df GVIF^(1/(2*Df))
## Species 482.320105 2 4.686341
## DHW4 3.455627 1 1.858932
## recent.meanDHW 2.111638 1 1.453148
## returnDHW6 1.734821 1 1.317126
## Species:DHW4 517.710799 2 4.770038In addition there is no colinearity between these core variables, which is good.
dredge.mod2 <- dredge(lm(ED50~Species*DHW4 +., data = dredge.dat1), rank = "AIC", m.lim = c(1,4))## Fixed term is "(Intercept)"#just look at the top scoring models
head(dredge.mod2, 10)## Global model call: lm(formula = ED50 ~ Species * DHW4 + ., data = dredge.dat1)
## ---
## Model selection table
## (Int) DHW rcn.mDH rDH Spc vSS DHW:Spc df logLik AIC delta weight
## 46 4.321 0.30680 0.08390 + + 8 1.479 13.0 0.00 0.391
## 44 6.986 0.21650 -0.2929 + + 8 1.129 13.7 0.70 0.275
## 16 6.607 0.12220 -0.1630 0.05856 + 7 -0.548 15.1 2.05 0.140
## 58 1.724 0.22950 + 1.4590 + 8 -0.365 16.7 3.69 0.062
## 14 5.212 0.17540 0.08732 + 6 -2.592 17.2 4.14 0.049
## 30 4.335 0.14940 0.06693 + 0.4729 7 -2.277 18.6 5.51 0.025
## 12 7.998 0.06906 -0.2910 + 6 -3.664 19.3 6.29 0.017
## 28 5.601 0.06648 -0.1862 + 0.7482 7 -2.720 19.4 6.40 0.016
## 11 8.734 -0.3401 + 5 -4.907 19.8 6.77 0.013
## 27 6.183 -0.2278 + 0.7871 6 -3.967 19.9 6.89 0.012
## Models ranked by AIC(x)#extracting the best fit model using AIC score
bestmodel <- get.models(dredge.mod2, 1)[[1]]
lm.dredge2 <- lm(bestmodel, data = dredge.dat1)
summary(lm.dredge2)##
## Call:
## lm(formula = bestmodel, data = dredge.dat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.46658 -0.17588 0.03949 0.13706 0.58307
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.32105 0.46482 9.296 7.51e-08 ***
## DHW4 0.30680 0.06293 4.876 0.000168 ***
## returnDHW6 0.08390 0.01889 4.442 0.000410 ***
## SpeciesP. meandrina 1.90207 0.61922 3.072 0.007300 **
## SpeciesP. verrucosa 1.94108 0.67319 2.883 0.010806 *
## DHW4:SpeciesP. meandrina -0.21359 0.08789 -2.430 0.027230 *
## DHW4:SpeciesP. verrucosa -0.18608 0.09230 -2.016 0.060924 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.272 on 16 degrees of freedom
## Multiple R-squared: 0.8185, Adjusted R-squared: 0.7504
## F-statistic: 12.03 on 6 and 16 DF, p-value: 3.728e-05ggeffects::ggeffect(lm.dredge2) %>% plot## $DHW4
##
## $returnDHW6
##
## $Species
When we now include the interaction recent.meanDHW drops out from the
best model. However, it is still important and results in a better
model.
lm.2 = lm(ED50~Species * DHW4 + recent.meanDHW + returnDHW6, data = dredge.dat1)
lm.5 = lm(ED50~Species * DHW4 + returnDHW6, data = dredge.dat1)
AICc(lm.2, lm.5)## df AICc
## lm.2 9 21.94864
## lm.5 8 23.32708anova(lm.2, lm.5)## Analysis of Variance Table
##
## Model 1: ED50 ~ Species * DHW4 + recent.meanDHW + returnDHW6
## Model 2: ED50 ~ Species * DHW4 + returnDHW6
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 15 0.87572
## 2 16 1.18410 -1 -0.30838 5.2821 0.03634 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Multiple interactions
lm.2 = lm(ED50~Species * DHW4 + recent.meanDHW + returnDHW6, data = dredge.dat1)
lm.6 = lm(ED50~Species * DHW4 * recent.meanDHW + returnDHW6, data = dredge.dat1)
lm.7 = lm(ED50~Species * DHW4 * returnDHW6 + recent.meanDHW, data = dredge.dat1)
lm.8 = lm(ED50~Species * returnDHW6 + DHW4 + recent.meanDHW, data = dredge.dat1)
lm.9 = lm(ED50~Species * recent.meanDHW + DHW4 + returnDHW6, data = dredge.dat1)
AICc(lm.2, lm.6, lm.7, lm.8, lm.9)## df AICc
## lm.2 9 21.94864
## lm.6 14 66.28199
## lm.7 14 61.85056
## lm.8 9 28.24851
## lm.9 9 29.57824None of the other interactions improved the model.
Model validation
library(performance)
check_model(lm.2)## Loading required namespace: qqplotr## For confidence bands, please install `qqplotr`.## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.## `geom_smooth()` using formula 'y ~ x'## `geom_smooth()` using formula 'y ~ x'## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 22 rows containing missing values (geom_text_repel).
r.squaredGLMM(lm.2)## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.## R2m R2c
## [1,] 0.8147412 0.8147412summary(lm.2)##
## Call:
## lm(formula = ED50 ~ Species * DHW4 + recent.meanDHW + returnDHW6,
## data = dredge.dat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38367 -0.08972 0.03304 0.09462 0.43907
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.78952 0.76071 7.611 1.58e-06 ***
## SpeciesP. meandrina 1.84012 0.55064 3.342 0.004460 **
## SpeciesP. verrucosa 2.19242 0.60783 3.607 0.002589 **
## DHW4 0.25490 0.06028 4.229 0.000729 ***
## recent.meanDHW -0.17568 0.07644 -2.298 0.036344 *
## returnDHW6 0.05255 0.02162 2.431 0.028076 *
## SpeciesP. meandrina:DHW4 -0.20679 0.07812 -2.647 0.018295 *
## SpeciesP. verrucosa:DHW4 -0.21270 0.08280 -2.569 0.021379 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2416 on 15 degrees of freedom
## Multiple R-squared: 0.8658, Adjusted R-squared: 0.8031
## F-statistic: 13.82 on 7 and 15 DF, p-value: 1.633e-05#Slopes comparison:
emtrends(lm.2, pairwise~Species, var = "DHW4")## $emtrends
## Species DHW4.trend SE df lower.CL upper.CL
## A. cf humilis 0.2549 0.0603 15 0.1264 0.383
## P. meandrina 0.0481 0.0581 15 -0.0758 0.172
## P. verrucosa 0.0422 0.0693 15 -0.1054 0.190
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## A. cf humilis - P. meandrina 0.20679 0.0781 15 2.647 0.0455
## A. cf humilis - P. verrucosa 0.21270 0.0828 15 2.569 0.0528
## P. meandrina - P. verrucosa 0.00591 0.0821 15 0.072 0.9971
##
## P value adjustment: tukey method for comparing a family of 3 estimateslibrary(jtools)
library(interactions)##
## Attaching package: 'interactions'## The following objects are masked from 'package:jtools':
##
## cat_plot, interact_plot, johnson_neyman, probe_interaction,
## sim_slopessim_slopes(lm.2, pred = DHW4, modx = Species, johnson_neyman = FALSE)## SIMPLE SLOPES ANALYSIS
##
## Slope of DHW4 when Species = P. verrucosa:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.04 0.07 0.61 0.55
##
## Slope of DHW4 when Species = P. meandrina:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.05 0.06 0.83 0.42
##
## Slope of DHW4 when Species = A. cf humilis:
##
## Est. S.E. t val. p
## ------ ------ -------- ------
## 0.25 0.06 4.23 0.00emtrends(lm.2, ~Species, var = "recent.meanDHW")## Species recent.meanDHW.trend SE df lower.CL upper.CL
## A. cf humilis -0.176 0.0764 15 -0.339 -0.0128
## P. meandrina -0.176 0.0764 15 -0.339 -0.0128
## P. verrucosa -0.176 0.0764 15 -0.339 -0.0128
##
## Confidence level used: 0.95sim_slopes(lm.2, pred = recent.meanDHW, modx = Species, johnson_neyman = FALSE)## Warning: recent.meanDHW and Species are not included in an interaction with one
## another in the model.## SIMPLE SLOPES ANALYSIS
##
## Slope of recent.meanDHW when Species = P. verrucosa:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -0.18 0.08 -2.30 0.04
##
## Slope of recent.meanDHW when Species = P. meandrina:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -0.18 0.08 -2.30 0.04
##
## Slope of recent.meanDHW when Species = A. cf humilis:
##
## Est. S.E. t val. p
## ------- ------ -------- ------
## -0.18 0.08 -2.30 0.04###Relative Importance of Predictors
library('relaimpo')
#calc.relimp(lm.2, type = "lmg")
calc.relimp(lm.2, type = "lmg")@lmg ## Species Species:DHW4 DHW4 recent.meanDHW returnDHW6
## 0.28829712 0.08514947 0.16291827 0.18219391 0.14721637sum(calc.relimp(lm.2, type = "lmg")@lmg)## [1] 0.8657751calc.relimp(lm.2, type = "lmg")@lmg / sum(calc.relimp(lm.2, type = "lmg")@lmg)## Species Species:DHW4 DHW4 recent.meanDHW returnDHW6
## 0.33299307 0.09835056 0.18817619 0.21044022 0.17003996For the output, we see that total proportion of variance explained by model is 86.6%. Species contributed 28.8% of variance, interaction of species * DHW4 contributed 8.5%, DHW4 contributed 16.3%, recent mean DHW contributed 18.2% and returnDHW6 accounted for 14.7% (total of these predictors = 86.6%, which is the total variance explained by the model)
Model Predictions
Fig3a = ggplot(mod.sep, aes(y = emmean, x = Reef, ymin = lower.CL, ymax = upper.CL, col = Species, shape = Species)) +
geom_pointrange(position = position_dodge(width = .4), size = .2) +
scale_y_continuous(expand = c(0,0), lim = c(6.1,8.6), breaks = c(6.5, 7, 7.5, 8, 8.5)) +
scale_colour_manual(values=col) +
labs(x="",
y="ED50 (95% CI)") +
theme_classic() +
theme(axis.line = element_blank(),
panel.border = element_rect(size = .5, fill = "transparent"),
legend.title = element_blank(),
legend.background = element_rect(fill = "transparent"),
legend.position = "bottom",
legend.direction = "vertical",
legend.text = element_text(face = "italic", size = 8, family = "Helvetica"),
legend.key.size = unit(3, units = "mm"),
legend.margin = margin(0,0,0,0, unit = "mm"),
axis.text = element_text(size = 8, family = "Helvetica"),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.title = element_text(size = 9, family = "Helvetica"),)
Fig3a
ggsave("2_figures/Fig3_ED50.pdf", plot = Fig3a, width = 70, height = 78, units = "mm")#DHW4
grid.DHW4 = with(dredge.dat1, list(DHW4 = seq(min(DHW4)-.1, max(DHW4)+.1, len = 50)))
em.DHW4 = emmeans(lm.2, ~DHW4|Species, at = grid.DHW4, type = "response") %>%
as.data.frame() %>%
mutate(var = "DHW4") %>%
rename(value = DHW4)
#returnDHW6
grid.returnDHW6 = with(dredge.dat1, list(returnDHW6 = seq(min(returnDHW6)-.1, max(returnDHW6)+.1, len = 50)))
em.returnDHW6 = emmeans(lm.2, ~returnDHW6|Species, at = grid.returnDHW6, type = "response") %>%
as.data.frame() %>%
mutate(var = "returnDHW6") %>%
rename(value = returnDHW6)
#recent.meanDHW
grid.recent.meanDHW = with(dredge.dat1, list(recent.meanDHW = seq(min(recent.meanDHW)-.1,
max(recent.meanDHW)+.1, len = 50)))
em.recent.meanDHW = emmeans(lm.2, ~recent.meanDHW|Species, at = grid.recent.meanDHW, type = "response") %>%
as.data.frame() %>%
mutate(var = "recent.meanDHW") %>%
rename(value = recent.meanDHW)
library(wesanderson)
col = wes_palette("GrandBudapest1", 4, type = "continuous")
DHW4.lm.plot = ggplot(em.DHW4, aes(x = value, y = emmean, fill = Species)) +
geom_point(data = dredge.dat1, aes(y = ED50, x = DHW4, col = Species, shape = Species), alpha = .8) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL), alpha = .2) +
geom_line(aes(col = Species)) +
scale_color_manual(values = col) +
scale_fill_manual(values = col) +
scale_y_continuous(expand = c(0,0), lim = c(6.2,8.6), breaks = c(6.5, 7, 7.5, 8, 8.5)) +
scale_x_continuous(expand = c(0,0), lim = c(5,11), breaks = c(6, 8, 10)) +
coord_cartesian(xlim = c(5.6,10.6), ylim = c(6.1,8.6)) +
labs(y = "ED-50 (95% CI)", x = "# events > 4°C-weeks") +
theme_classic() +
theme(axis.line = element_blank(),
panel.border = element_rect(size = .5, fill = "transparent"),
legend.position = "none",
axis.text = element_text(size = 8, family = "Helvetica"),
axis.title = element_text(size = 9, family = "Helvetica"),
axis.title.y = element_blank())
returnDHW6.lm.plot = ggplot(em.returnDHW6, aes(x = value, y = emmean, fill = Species)) +
geom_point(data = dredge.dat1, aes(y = ED50, x = returnDHW6, col = Species, shape = Species), alpha = .8) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL), alpha = .2) +
geom_line(aes(col = Species)) +
scale_color_manual(values = col) +
scale_fill_manual(values = col) +
scale_y_continuous(expand = c(0,0), breaks = c(6.5, 7, 7.5, 8, 8.5)) + #lim = c(6.2,8.6)
scale_x_continuous(expand = c(0,0), breaks = c(3,6, 9, 12)) +
coord_cartesian(xlim = c(1.85, 13.09), ylim = c(6.1,8.6)) +
labs(y = "ED-50 (95% CI)", x = "Return time events > 6°C-weeks") +
theme_classic() +
theme(axis.line = element_blank(),
panel.border = element_rect(size = .5, fill = "transparent"),
legend.position = "none",
axis.text = element_text(size = 8, family = "Helvetica"),
axis.title = element_text(size = 9, family = "Helvetica"),
axis.title.y = element_blank())
recent.meanDHW.lm.plot = ggplot(em.recent.meanDHW, aes(x = value, y = emmean, fill = Species)) +
geom_point(data = dredge.dat1, aes(y = ED50, x = recent.meanDHW, col = Species, shape = Species), alpha = .8) +
geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL), alpha = .2) +
geom_line(aes(col = Species)) +
scale_color_manual(values = col) +
scale_fill_manual(values = col) +
scale_y_continuous(expand = c(0,0), breaks = c(6.5, 7, 7.5, 8, 8.5)) + #lim = c(6.2,8.6)
scale_x_continuous(expand = c(0,0), breaks = c(4,5,6)) + #lim = c(27.2,28.8),
coord_cartesian(xlim = c(3.44,6.6), ylim = c(6.1,8.6)) +
labs(y = "ED-50 (95% CI)", x = "mean maxDHW 2016-20") +
theme_classic() +
theme(axis.line = element_blank(),
panel.border = element_rect(size = .5, fill = "transparent"),
legend.position = "none",
axis.text = element_text(size = 8, family = "Helvetica"),
axis.title = element_text(size = 9, family = "Helvetica"),
axis.title.y = element_blank()) gpt = theme(axis.title = element_blank(),
plot.margin = unit(c(1,1,1,1), "mm"))
ggarrange(Fig3a + theme(legend.position = "none"),
DHW4.lm.plot + gpt,
returnDHW6.lm.plot + gpt,
recent.meanDHW.lm.plot + gpt,
ncol = 4, align = "hv")
ggsave("2_figures/Fig3_model_v3.pdf", width = 180, height = 53, units = "mm")